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I. INTRODUCTION AND SUMMARY 



When correlations and fluctuations become important in the dynamical evolution of a 
many-body system and mean field approximations are not sufficient, molecular dynam- 
ics methods are invoked frequently. Molecular dynamics means that the constituents 
(molecules) of the many-body system are represented by few classical degrees of freedom 
(center of mass position and momentum, angles of rotation etc.) and interact through po- 
tentials. The equations of motion (Newton's or Hamilton's type) are solved numerically. 
The interactions which are used range from purely phenomenological to sophisticated ab 
initio quantal potentials. The major advantage of molecular dynamics simulations is that 
they do not rely on quasi-particle approximations but include both, mean-field effects and 
many-body correlations. Therefore they can provide insight into complex systems with 
correlations on different scales. 

During the past decade one can observe an increasing interest in the dynamics of many- 
fermion systems in which correlations are important. In nuclear as well as in atomic physics 
collisions of composite fermion systems like nuclei or atomic clusters are demanding many- 
body models which can account for a large variety of phenomena. Depending on the energy 
one observes fusion, dissipative reactions, fragmentation and even multi-fragmentation, va- 
porisation or evaporation and ionization. Also phase transitions in small systems are of 
current interest. 

Classical molecular dynamics is applicable if the de Broglie wave-length of the molecules 
is small compared to the length scale of typical variations of the interaction, otherwise the 
quantal uncertainty relation becomes important. If the molecules are identical fermions or 
bosons the de Broglie wave-length should also be small compared to the mean inter-particle 
distance in phase space, otherwise the Pauli or Bose principle is violated and the model will 
have wrong statistical properties. 

For nucleons in a nucleus, for example, both conditions neccessary for classical mechanics 
are not fulfilled. The same holds for electrons in bound states or at high densities and low 
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temperatures. Nevertheless one would like to utilize the merits of a molecular dynamics 
model for indistingushable particles in the quantum regime. 

This article reviews attempts that are made to combine Fermi-Dirac statistics with a 
semi-quantal trajectory picture from the viewpoint of the quantal time-dependent varia- 
tional principle. The closest quantum analogue to a point in single-particle phase space 
representing a classical particle is a wave packet well localized in phase space. The analogue 
to a point in many-body phase space representing several classical particles is a many-body 
state which is a product of localized single-particle packets. If the particles are identical 
fermions this product state has to be antisymmetrized, in the case of bosons it has to be 
symmetrized. 

For models which are formulated in terms of trial states and a Hamilton operator, both, 
static and dynamical properties can be obtained from appropriate quantum variational prin- 
ciples. Ground states can be determined with the help of the Ritz variational principle and 
the equations of motion can be accessed through the time-dependent variational princi- 
ple which allows to derive approximations to the time- dependent Schrodinger equation on 
different levels of accuracy. 

In chapter |I| we first discuss the time-dependent variational principle in general and then 
show how it works for various trial states and how classical mechanics can be obtained from 
quantum mechanics by an appropriate choice of the dynamical variables. 

Using wave packets automatically guarantees that the Heisenberg uncertainty principle 
is not violated by the model. That is actually a great problem in classical simulations of 
Coulomb systems, where for instance the hydrogen atom for a point like proton and electron 
is infinitely bound. 

Using antisymmetrized many-body trial states automatically guarantees that the Pauli 
exclusion principle is respected by the model. Phenomena like shell structure or Fermi-Dirac 
statistics emerge in a natural way and are discussed in chapter and |V|. 

Although antisymmetrized product states of single-particle gaussian wave packets possess 
already major relevant degrees of freedom, some processes like disintegration of wave packets, 
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which can occur in coordinate space (evaporation, capture or tunnehng) and in momentum 
space (large momentum transfer due to colhsions), are poorly described by the equations of 
motion. One therefore often represents the system with a mixture of trial states between 
which random transitions may occur. This branching procedure, which is employed in atomic 
as well as in nuclear physics, will be explained in Sec. [11 C 2\ 

In chapter |I^ models used in nuclear and atomic physics are reviewed from the general 
point of view of chapter J| and |T|. When energetic collisions between heavy atomic nu- 



clei became available classical molecular dynamics models were developed to describe the 
various phenomena observed. To simulate the effect of antisymmetrization on the classical 
trajectories many authors add to the hamiltonian a two-body "Pauli potential" which is 
supposed to keep fermions apart from each other in phase space. 

Quantum Molecular Dynamics (QMD) attributes to each fermion a gaussian wave packet 
with fixed width instead of a point, but still uses a simple product state for the many-body 
wave function and therefore obtains classical equations of motion with two-body forces acting 
on the centroids of the wave packets. These Newtonian forces are supplemented by random 
forces which simulate hard collisions. In these collision terms Pauli blocking is included. 
The statistical properties are however mainly those of distinguishable particles. 

Antisymmetrized Molecular Dynamics (AMD) also uses gaussian wave packets with fixed 
width, but antisymmetrizes the many-body state. Like in QMD a collision term is added 
to account in a phenomenological way for branching into other Slater determinants. In 
Fermionic Molecular Dynamics (FMD) in addition the width degree of freedom is considered. 
This non-classical degree of freedom is important for phenomena like evaporation and it also 
plays an important role in statistical properties of the model. 

Section [IV C| explains the different use of the term Quantum Molecular Dynamics in 
atomic and nuclear physics. It then refers to applications of trajectory calculations for 
individual electrons and ions where the density of the electrons is too large to neglect their 
fermion character. Quantum branching in the atomic context is also briefly discussed. 

Although designed for non-equilibrium simulations like collisions, molecular dynamics 
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models are also used to simulate systems in thermal equilibrium. In chapter their sta- 
tistical properties are investigated by means of time averaging. As applications the nuclear 
liquid-gas phase transition or the hydrogen plasma under extreme conditions are discussed. 

In conclusion, the review shows that it is possible to extend the classical trajectory picture 
to identical fermions by means of localized wave packets. When the phase space density 
increases the classical notion of positions and momenta pk as mean values of narrow wave 
packets has to be reinterpreted as parameters which identify an antisymmetrized many- 
body state \Q) = I ri,pi, r2,p2, ■ ■ ■ )• When the individual packets overlap in phase space 
fk and Pk can no longer be identified with the classical variables. Calculating all observables 
with the many-body state \Q) as quantum expectation values, e.g. 'H{ri,pi,r2,P2, ■ ■ ■) = 
{Q \ H \ Q) / {Q \ Q) , and not misinterpreting them as classical expressions, includes the 
Pauli exclusion principle and Heisenberg's uncertainty principle in a natural and correct 
way. In the equations of motion the exclusion principle causes a complicated metric in the 
A^-particle parameter space in the sense that canonical pairs of variables can only be defined 
locally. 

The antisymmetrization of localized wave packets, which brings together the Pauli prin- 
ciple, Heisenberg's uncertainty relation and the classical trajectory picture, leads to many, 
some times unexpected, quantal features. Nevertheless one has to be aware that one is still 
dealing with a very simplified trial state and other degrees of freedom may be important. 
Especially when the interaction is not smooth across a wave packet it may want to change its 
shape to one which is not in the allowed set, for example split with certain amplitudes into 
different parts which after some time evolve independently. Or more generally, branching 
into other trial states away from the one which follows the approximate time evolution of a 
pure state leads to a mixture of antisymmetrized wave packets. The consistent treatment of 
this aspect of quantum branching needs further attention in the literature. 

Common to all models discussed is the anticorrelation between its degree of consistent 
derivation and its computational effort. 
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II. TIME-DEPENDENT VARIATIONAL PRINCIPLE 



A. General remarks 

The time evolution of a state in quantum mechanics is given by the time-dependent 
Schrodinger equation [] 

^J^\m) = H\m) , (1) 

where H is the hamiltonian and | ) the many-body state which describes the physical 
system. This equation of motion can be obtained from the variation of the action ( |Kerman| 
and Koonin| , |1976| ; [Kramer and Saraceno| , [198 1| ; prozdz et al.\ , [I986| ; [Broeckhove et al] , |1988| ) 



A'= r'rft/:'(v^(t)*,vl>(t),vi>(t)) , (2) 

Jtl 

keeping the variations fixed at the end points: S'^iti) = 5\l/(t2) = S'^*{ti) = (5\l>*(t2) = 0. 
The Lagrange function 

c'mtr, m, ^ w) = ( m i i ^ w ) - ( m i h \ m ) (3) 

is a function of the dynamical variables, denoted by the set Usually they are chosen to 

be complex so that and may be regarded as independent variables. If the mapping 
of on the many-body state | ) is analytic, | ) depends only on and the 
hermitian adjoint state (\&(t)| only on the complex conjugate set '$*{t). The Lagrange 
function depends on the first time derivatives ^(t) through ^ | \l/(t) ). 

In general, the set \&(t) contains infinitely many dynamical degrees of freedom, for exam- 
ple the complex coefficients of an orthonormal basis or the values of a A^-body wave function 
on a 3A^-dimensional grid in coordinate space. 



^ Throughout the article all operators are underlined with a tilde, e.g. H, and expectation values 
are denoted by caligraphic letters, e.g. 7i. If not needed explicitly h is taken to be one. 
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The Lagrange function is appropriate if the set of variables contains a 

complex overall factor to | ^(t) ) which takes care of norm and phase. If | ) cannot be 
normalized by means of the set \l'(t) the Lagrange function defined in Eq. (0) should be 
used. We shall employ both forms. 

For linearly independent variables ^(t) and \E'*(t) the variation of \E'*(t) yields for the 
extremal action 

= 5A'=rdt |z^ . (4) 

Jti dt ~ 

If I \&(t) ) represents the most general state in Hilbert space, the variation ( 5\I'(t) | is unre- 
stricted and Eq. @ can only be fulfilled if 

0|-^)l*W) = O, (5) 

which is just the Schrodinger equation (|1|). 
The variation of \I'(if:) yields 

= SA' = iyt {^^Iit)\^J^-H\S^it)) , (6) 

which after partial integration over t with fixed endpoints | S'i/{ti) ) = | 5\l/(t2) ) = results 
in the hermitian adjoint of the Schrodiger equation 

(*Wl(-4-?)=0. (7) 

As is linearly independent of \l/*(t), both have to be varied. Therefore the time 
integral which was missing in early attempts, used by ( [Frenkel| , |1934] ) and others, is necessary. 



The reason for the reformulation of a differential equation as a variational principle is 
of course the anticipation that a suitably chosen restriction in the dynamical variables will 
lead to a useful approximation of the full problem. For more general considerations on 
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Balian and 



Veneroni], |1988|) and references therein. 
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1. Restricted set of variables 



Let a restricted choice of variables be denoted by the complex set Q'{t) = 
{qo(t) , qi{t) , q2{t) , ■ ■ ■} which specifies the many-body state \Q'{t)). It is presumed that 
qo{t) is always a complex overall factor in the sense 

I Q'it) ) = qoit) I Q{t) ) = qo{t) I q2{t), ■■■) . (8) 

Furthermore, the time-dependence of | Q'{t) ) is supposed to be only implicit through the 
variables Q'(t). 

It should be noted that the manifold | Q'{t) ) is in general only a subset of the Hilbert 
space and needs not to form a subspace. 

Variation of the action (|^) with the Lagrange function (|^) with respect to Q'* [t) leads to 



= 6A'= f'dtY^dqlit) f/^ 



As Sql(t) are arbitrary functions the action is extremal if the following equations of motion 
are fulfilled 

^ ^ ■{Q'it)]] (^^-H] \Q'{t)) = (10) 



dqt J \ dt 



or 



^j:Cl,q,=^{Q'{t)\H\Q'{t)) , (11) 

A* qu 

where 

Different from the unrestricted variation one obtains equations of motion for the complex 
parameters which in turn define the time evolution of the trial state | Q'{t) ) in Hilbert space. 

After partial integration over time variation with respect to Q'{t), with fixed end points, 
^qu{ti) = 5qy{t2) = 0, results in equations of motion which are the complex conjugate of 
Eqs. (H - (0). 
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In the following, for sake of simplicity, the explicit indication of the time dependence is 
sometimes omitted. 

The time evolution of go can be expressed in terms of the other variables. For that the 
Lagrangian CJ may be written as 

at 

= ^gS9o(QIQ) + go*9o(QIQ)/:(Q,Q*,Q,Q*)-^^(go?o(g|g)) , (13) 

where the set Q does not contain go anymore (see Eq. (||)) and a new Lagrangian 
L{Q\Q,Q\Q) is defined by 



C{Q\Q,Q\Q) = l- 



2V {Q\Q) J {Q\Q) 

CoiQ*,Q,Q*,Q)-niQ*,Q) (14) 



with 



iQ>-li<3)=i:4.^io> -d (15) 

The new Lagrange function C contains the norm explicitly and is made real by subtracting 
the total time derivative. T-C{Q*,Q) is the expectation value of the hamiltonian H and will 
be called Hamilton function. 

It is easy to verify that the solution of the equation of motion ([To|) for go is 

go(t) = ^exp h f'dt'C{t')} (16) 

and analogue for the complex conjugate go(i)- Thus, the variational freedom of an overall 
factor go(t) is used by the time-dependent variational principle to provide a state | Q'{t) ) 
with a time-independent norm and an the additional phase /* dt'C{t') {C is real by construc- 



tion). Furthermore, insertion of the solution (|l^) into C shows that along the trajectory 



/:'(g'*(t),g'(t),g'(t)) = o, (i7) 

irrespective of the choice of the remaining degrees of freedom in \Q). 



11 



The equations of motion ([T0|) and their complex conjugate ones can of course also be 
expressed as Euler Lagrange equations 

ddC _^ ddC _dC' _^ 

dt dq* dq* dtdq^dqu 

For the remaining variables {g*, q^] u ^ 0} they can be written in terms of the new Lagrange 
function £ as 

d dC dC ( d dC dC\ , 

qMQ\Q) + 



dtdq* dq* \dtdq* dq* 

j^i^M Q\Q))^^- i^<lom + q*oqoC)^^§^ = (19) 

and the analogue complex conjugate ones. The last two terms vanish when the general 
solution (|l^) for go is inserted. Hence for u ^ the Euler Lagrange equations with C as 
given in ( p!4D are equivalent to the ones with C. 
Therefore, the action 

A= r dt C{Q*,Q,Q*,Q) , (20) 
which one also finds often as a starting point ( [Kerman and Koonin , 1976 ; Kramer and 



Saracenc] , 1981 ; Urozdz et al. , 1986 ; Broeckhove et a/] , |1988| ), is the appropriate one if the 



trial state | Q ) is not normalized and the phase is disregarded]^. This form will turn out to 
be more convenient when dealing with antisymmetrized states in section p| . 

One should however keep in mind that the simple inclusion of go provides automatically 
norm and phase. 

The Euler Lagrange equations, which result from variation of the action (^) 

d dC d C ^ d dC dC 

dtdql dq* dtdq^ dqu 



^£ may also be written with normalized trial states \ Q) / \/ {Q\Q) ■ The only difference is a total 
time derivative emerging from so that the Lagrange functions are equivalent. 
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can be written in terms of the Hamilton function Ti = {Q\H\Q)/{Q\Q) as generalized 
Hamilton's equations 

«E'^M-'?- = S and -^Y.^lu€ = ^ ■ (22) 
The nonnegative hermitian matrix C depends in general on Q and is given by 

(010) {Q\Q) (010) ^ ' 

It plays the role of a metric on the manifold of dynamical variables Q. 

There is no need to assume that the set Q contains only complex variables. But as will 
be seen later in the context of classical mechanics, the real and imaginary parts of play 
the role of canonical pairs of variables. On the other hand for quantum mechanics it is only 
natural to work with complex variables as the quantum state \Q) is necessarily complex. 
The assumption of an analytic mapping of Q onto the trial state | Q ) is also not compulsory 
but makes the equations usually more transparent. 



2. Deviation from the exact solution 

The static analogue to the time-dependent variational principle is the Rayleigh-Ritz 
variational principle for the energy. Here one can show that the overlap between the true 
ground state of the hamiltonian and the trial state \Q') is increasing when the energy 
{Q' \ H \ Q') is decreasing. Since any additional degree of freedom lowers the energy, more 
degrees of freedom imply always an improved description. General statements like that 
cannot be made for the time-dependent variational principle, although one would believe 
that more degrees of freedom will result in a time evolution which is closer to that of the 
Schrodinger equation. 

In order to understand better in which sense the time-dependent variational principle 
optimizes the evolution of the trial state, we consider the deviation of the approximate 
solution from the exact one, which develops during a short time 6t 
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FIG. 1. Sketch of the manifold of trial states with the approximate time evolution 
(solid line) and the solution of the Schrodinger equation (dashed line). The error 
I A(t, 6t) ) is orthogonal to all tangent states Q' ). 

I A(t, 6t))= I ^exarf(t + St))- \ Q\t + 5t)) = exp{-tH6t) \ Q\t) ) - \ Q\t + 5t) ) 

= -HH-iY^ qu^] I Q'{t) )5t + order(5t2) . (24) 

The equations of motion (0) which result from the time-dependent variational principle 
demand that 

/ d 



dq* 



{Q\t)\ \A{t,6t)) = 0, (25) 



which means that the deviation is orthogonal to all tangent states ^-| Q' ) , see Fig. |l|. In 
other words the approximate equations of motions evolve Q'{t) to that point Q'{t + 6t) 
in the manifold where any small change in all possible directions increases the distance 
( A(t, 6t) I A(t, 6t) ) between true | "^exactit + 6t) ) and approximate solution | Q'(t + 6t)). 
The minimum condition 

0-^^^(A(Mt)|A(Mt)) 

^ ^ {Q'it + 6t)\] \A{t,6t)) 



dq*{t + 5t) 
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^ ( Q'{t) I ) I A(t, 6t) ) + order((5t^) . (26) 



dqtit) 



is fulfilled because of Eq. (25). For variations with respect to the complex conjugate of 
Eqs. (p5|) and (p6|) have to be used. Of course after many time steps 6t the deviation may 
become large. 

3. Poisson brackets and canonical variables 

With the help of C it is possible to introduce Poisson brackets. Using the equations of 
motion (^) the time derivative of the expectation value of a time-independent observable 
B can be written as 

d {Q\B\Q) f.^dB . dB 



dt (QIQ) ^Vdql ^"dq, 
i 



^^dq, ^'"dqt dq, d q*. 



= :{H,B}. (27) 
The real and imaginary parts of q^ form pairs of canonical variables if C is the unit matrix. 



In the following section |II B| two examples are given for this case, the Schrodinger equation 
and Hamilton's equation of motion. 

In the general case, in which C^^, is not diagonal and depends on Q, pairs of canoni- 
cal variables exist locally according to Darboux's theorem ( Arnol'd , 1989| ). One possible 
transformation is given by 

dr, = E C;/' (Q* , Q) dg. and dr^ = ^ ^^l (Q* , Q) ■ (28) 

Written with the new variables and r* the Poisson bracket (^) takes the form 

dU dB dB OH 



so that {rij,,r*^) form pairs of canonical variables. The problem is, however, that in 
nontrivial cases the transformation ( p8|) cannot be written in a global way as R{Q) = 
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{f^oi^O: Q'l! ■ ■ ■'^1(^0, Q'l, ■ ■ ■ ■ ■} the Hamilton function H or the observable B cannot 
be expressed in the new variables 

n{R\R)=n{Q\Q) , B{R\R) = B{Q\Q) . (30) 

A set of canonical pairs which are real is given by 

i=(r; + r,) and vr, = 



Pi^ = ^ + '^m) and 71-^ = ^ - r^) , (31) 



which yields the standard Poisson brackets 



{„,e}.5:f|^|i-|i|^), (32) 



where Ti and B are now functions of and vr^. Besides the trivial examples discussed in 
Sec. [II B| we give one nontrivial example in Sec. [Ill A 5 . 



The mere fact that according to Darboux's theorem canonical pairs (p^, vr^) exist, ad- 
mits to guess or make an ansatz for the Hamilton function 7Y(pi, • • • ? tti, 7r2, . . .) and use 
Hamilton's equations of motion. But with such a guess the connection to the trial state 
I Q ) is lost and the physical meaning of p^ and tt^ is obscured. This will become obvious 
when we discuss trial states for indistinguishable particles in Sec. PTT] . 

4- Conservation laws 



An expectation value is conserved if ( Broeckhove et al. , 1989| ) 



^V9g^, dq* dq^ ^ dq*J 
Hence the energy Ti itself is always conserved by the equations of motion, provided they are 
derived from the variational principle. This is completely independent on the choice of the 
trial state. 

In the following it is shown how to identify other constants of motion and how the trial 
state has to be chosen in order to ensure desired conservation laws. For that we consider a 
unitary transformation with the hermitian generator G 
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U = exp (^ie , e real . (34) 

If U maps the set of trial states onto itself 

U\Q') G HQ')}, (35) 

then the special infinitesimal variation | Q'(t) + 6Q'{t) ) = exp{i6e(t) G} \ Q'(t) ) of the ac- 
tion (H) yields 

= j^*' dt ( Q\t) I exp{-i 6s{t) G}(^i^- exp{t 6s{t) G} \ Q'{t) ) (36) 
= dt 5e{t) \j^{Q\t) I G I Q'it) ) - (Q'it) 1 1 [h,g] \ Q'it))\ 



Hi 

-|-total time derivative + order(5e^ 



As 6e{t) is arbitrary and vanishes at the end points one obtaines 

|a = ^(Q'wieiQ'w) 

= {W, 5} = ( Q'{t) I ! [h, g] I Q'(t) ) . (37) 

That means that for this class of generators the generalized Poisson bracket is just the 
expectation value of the commutator with iH. 

Relation (P?]) is very useful for two reasons. First, if G commutes with the hamiltonian 
H and exp{ieG) \Q'i) = \ Q2 ) then {Q'{t) \ G\ Q'{t) ) is automatically a constant of motion. 
Second, this relation is an important guidance for the choice of the trial state \Q')- If one 
wants the model to obey certain conservation laws then the set of trial states should be 
invariant under unitary transformations generated by the constants of motion. For example 
total momentum conservation implies that a translated trial state is again a valid trial state. 
Conservation of total spin J = L + S is guaranteed when a rotation of the trial state in 
coordinate and spin space results again in a trial state. 

Relation (|37|) also sheds some light on the quality of the variational principle. It says that 
under the premises that exp{isG) does not map out of the set of trial states the expectation 
value Q(t) of G develops for short times like the exact solution. From Eq. (|37| ) follows that 
along the trajectory | Q'{t) ) the time derivative of Q{t) equals the exact one. 
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B. From quantum to classical mechanics 



This section demonstrates that the time-dependent variational principle, discussed in 
general in the previous section, represents a method to go in a smooth way from quantum 
physics to classical physics by appropriately choosing the dynamical degrees of freedom in 
the trial state. In contrast to Ehrenfest's theorem this method also works for identical 
particles and in finite dimensional spin-spaces. 

1. Quantum mechanics 

As a first illustration let us represent the trial state in terms of an orthonormal basis 
I n ) in many-body space. We may write a general state \Q') as 

\Q') = |pi,P2, ■ ■ ■ ,vri,7r2, ■ ■ ■ ) = ^ 4^*^^" + \ ^) = Y1 (^n\n) , (38) 

where the complex amplitudes are written in terms of their real and imaginary parts p„ and 
Tin- It is easy to verify that the Lagrange function for normalizable states defined in Eq. 
is given by 

1 d i 

^' = J2 2 ~ ^"^"^ ^ rft ^ 4 ^ ~ /'s, ■ ■ ■ , TTi, 7r2, ■ ■ ■) • (39) 

The real Hamilton function, expressed with the real and imaginary parts of the matrix 
elements Hnk = {f^l H \ k) , is bilinear in p and tt 

H{pi,P2, ■ ■ ■ , VTi, 7r2, ■ ■ ■) = ^ XI liPkPn + 7rfc7r„) ReHkn + {-T^kPn " PkT^n) Imiffcn] • (40) 

fc,n 

The Euler Lagrange equations 

d dC dC , d dC dC 

dt O Pn O Pn dt OTTn O Tin 

yield 

d an d an 

-n^n = — and —pn = ^ — . (42) 

dt a Pn dt ann 



There is a very important message to be learned from this little exercise: Eqs. (^2]) look 
exactly like Hamilton's equations of motion in which (p^, 7r„) are pairs of canonical vari- 
ables and 7Y(pi, p2, - ■ ■ , vri7r2, ■ ■ ■) is the Hamilton function, bilinear in the coordinates and 
momenta of the system. But these seemingly classical equations are just a representation 
of the Schrodinger equation as can easily been seen by rewriting the two real Eqs. (E2p in 
terms of the complex coefficients c„ = -^{Pn + ^7r„) as one complex equation, namely 

«^c„ = = ^ Hnk Ck or i^\Q') = H\Q') . (43) 

n k 



The mere fact that the equations of motion (|4^) appear in a classical form does not neces- 
sarily imply that the system is classical and for example violates the uncertainty relation or, 
in case of indistinguishable fermions, Fermi-Dirac statistics. This will be discussed in detail 
in Sec. 0. 



The symplectic structure of Eqs. (|42|) is fundamental to all energy conserving dynamical 
theories, classical, quantum or quantum field theories, see e.g. ( |Katz| , |1965| ). Only the 
physical meaning of the dynamical variables p„ and 7r„ and of the Hamilton function Ti, 
determine which kind of physical system one is dealing with. 

In the example (^) p„ and 7r„ are not the 6A^ positions and momenta of the particles 
but infinitely many variables which specify the many-body state \Q')- A truncation to a 
finite number results in an approximation of the exact time evolution. The quality depends 
on how well the selected set of basis states, { | can represent the part of the Hilbert 
space which is occupied by the physical system under consideration. 

2. Classical mechanics 



A second example shows ( [Heller] , |1975| ) that one gets the classical Hamilton equations 



of motion for the 6A^ positions and momenta from the time-dependent variational principle 
by choosing the following dynamical variables. The normalized trial state | Q ) is set up to 
describe A^ distinguishable particles which are localized in phase space. 
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\Q) = 1^2,^2) \ rN,PN) (44) 

The single-particle states | , ) in the direct product are taken to be the closest quantum 
analogue to classical particles, namely gaussian wave packets of minimum uncertainty, i.e. 
coherent states ( [Klauder and Skagerstam] , |1985D : 



a^kfc.Pfc = exp<^ ^iVk^\ 45) 

\k\rk,Pk) = y—j exp<^ ^ irk{k-pi,)\. (46) 

The dynamical variables and are just the mean values of the position and momentum 
operator, respectively 

rk = {rk3k\^\rk3k) > Pk = {rk3k\i\rk3k) ■ (47) 

The width parameter qq is a real fixed number here. In later applications it will also be 
taken complex dynamical variable. 



The evaluation of the Lagrange function (|1^) is again simple and yields 



N 

C = -J2r,^k-n, (48) 

k=l 

where the Hamilton function Ti = {Q \ H \ Q) 

TV / -2 o \ N 

^ = E hr^ + I + E i^k,Pk\ ®{fi,Pi\V{l,2)\n,pk)® \ri,pi) (49) 

is the expectation value of the hamiltonian 

H = T.^+ E ^(^.0- (50) 

1=1 ^^l k<l=l 



The Euler Lagrange equations 



d dC dC ^ ^ dU 

d d C dC J, dH 

rk = -TT^ , (52) 



dtdpi, dpk dp, 

result in 

20 



dn , d^ dU 
—Pk = and — rfc = . (53) 

dt ovk dt apf. 

In a situation where classical mechanics holds, i.e. the wave packets are narrow enough so 
that one can approximate the expectation value 

{Q\H\Q) =niri,r2,---,Pi,P2,---) ^ H(fi,r2,---,Pi,P2,---) (54) 

by replacing the momentum and position operators in the hamiltonian by their respective 
mean values of the wave packets, Eqs. ( [53| ) become the classical Hamilton's equations of 
motion. 

Up to this point the time-dependent variational principle leads to the same results as 
Ehrenfest's theorem which usually establishes the connection between quantum and classical 



systems. As will be demonstrated in section |ITI| the prescription to replace the operators in 



the Heisenberg equation by their mean values does not work for indistinguishable particles. 
But the time-dependent variational principle for an antisymmetrized trial state | Q ) provides 
the molecular dynamics equations for identical fermions. For bosons one would of course 
use a symmetrized state. 

The canonical pair {r k,p of real position and momentum can be combined to a complex 
variable = y^^^fc + ^\J^Pk- The Hamilton equations (|53D written in terms of Zk and 
zl take the form 

d d 

i-^Zk = ^^(^1' 4r ■ ■ ,zi,Z2, - ■ ■) , (55) 

which formally looks like the Schrodinger equation (^). 

This and the previous example show that one cannot decide from the form of the equa- 
tions of motion alone whether the system, they are describing, is classical or quantal. Fur- 
thermore, the time-dependent variational principle can provide both, classical and quantal 
many-body equations of motion, depending how the trial state \Q) is chosen. In the fol- 
lowing third example an intermediate situation is sketched, where some quantum effects are 
included. 
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3. Semi-classical, semi-quantal 



A third example which has part of quantum effects included is a trial state in which 
the width parameter a is a dynamical variable and complex, a = + iaj ( [Isue anH 
Kujiwara , 1991 ). Take the A^-body trial state to be a product state 



\Q)= ki ) ® I g2 ) ® ■ ■ ■ ® I g^v ) (56) 

of single-particle states, which are gaussians 

(x = (x |n,pj,az) = 27r-— exp <^ VipiX+i<pi) (57) 

\ + ai J y lai j 

characterized by their mean positions r;, mean momenta p;, widths ai and phases 

For a single-particle hamiltonian which contains in addition to the kinetic energy a 
harmonic oscillator potential 

N N (k'^il) 1 \ 

^HO = E k^oil) = E ( ^ + -^muj^m] ■ (58) 
the equations of motion are 

d _^ Pl d _^ 2-* /r-r\\ 

-T:ri = — , -tPi = -1^1 ^ (59) 
dt mi at 

d i . 2 2 d pf 3 am mi 

ttO/ = tmiuai, —(j)i = -- - uj ri . (60) 

at mi at 2mi 2mi\ai\'^ 2 

These equations, although looking classically for fi and p;, represent the exact solution of 
the Schrodinger equation, provided the wave function is a gaussian wave packet at time zero. 
The centers of the wave packets as well as the mean momenta oscillate harmonically with 
the frequency uj. Due to the time dependence of the widths the packets also breath but with 
twice the frequency. The solution is fully quantum mechanical, although described by only 
a few (classically looking) parameters. Free motion without a potential [uj = 0) is of course 
also exact. For general potentials the trial state (|57| ) may serve as an approximation if 
locally, in the region where {x \qi) does not vanish, the potential is well represented by a 
harmonic oscillator. 
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In general one can say, that a trial state provides the exact solution of the Schrodinger 
equation if the action of the Hamilton operator H on the trial state can be expressed in 
terms of parameters and first derivatives with respect to the parameters like in the following 
example 



\ dr d(b. 



(61) 



and 



' d ^ d 
2« 7^ h 2zp ■ — - z 

. a aj or 



p + 3 



4 + 



d_ 
'dd 



q) 



(62) 



For the harmonic oscillator the Schrodinger equation takes the form (index / omitted) 



dqu 
' p\ d 
m j dr 



(63) 
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From Eq. ( |63D the equations of motion ( |59D and ( |60D follow at once. Because gaussian wave 
packets are an exact solution of the Schrodinger equation for these one-body hamiltonians, 
also the respective product states are an exact solution for the corresponding many-body 
problem. Moreover, since antisymmetrization and symmetrization commute with the exact 



time evolution, the equations of motion (pOD remain the same for antisymmetric product 
states (identical fermions) and symmetric product states (identical bosons). 

The time- dependent width parameter a, which describes the variances in coordinate 
and momentum space, provides the first nonclassical degree of freedom in the parameter 
manifold. It completes the classical equations (^) to the full quantum solution for spherical 
harmonic oscillator potentials as well as for free motion. In section ^ it is shown that also 
in thermodynamic considerations the inclusion of finite widths leads to quantum statistics. 
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C. Further remarks 



In the previous Sec. p.lB| several examples are presented to get acquainted with the 
numerous aspects of the time-dependent variational principle and its function as a bridge 
between quantum and classical physics. We add now further remarks which are again of 
general type and useful to understand the different models referred to later. 



1. Self consistency and nonlinearity 

The equations of motion (p!Tl) (or (|22|)) determine the time evolution of the parameters 
qu{t) in parameter space. But, even if these parameters have an intuitive physical meaning, 
they first of all determine the trial state | Q'{t) ) from which all physical observables have to 
be calculated in a quantum fashion. Therefore, we derive here the corresponding equation 
of motion for | Q'{t) ) in Hilbert space. 

The Hamilton operator Hq, which evolves the trial state | Q'{t) ) in time according to 
the equations of motion (|ll]), has to fulfill 



It is given by 



H, = t(\Q'){Q'\ - \ Q'){Q'\) + Uq')({Q'\Q')-{Q'\Q')){Q'\ , (65) 



2 



where | Q' ) stands for 



\Q') = i\Q') = -^i:^ I Q' ) C"' 1^ (66) 

dt j;idq^ ^ dq* 

and analogue for {Q' \ ■ Thus Hq itself depends on Q'* and Q' . Using the fact that | Q') 
is always normalized it is easy to show that defined in (|65|) is the generator of the 
approximate time evolution 

tj^\Q') = HoiQ'\Q') \Q') . (67) 
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The equation of motion ( |67D for the trial state in Hilbert space is the counterpart for the 
equation of motion (|Tl|) in parameter space 

= Y.C'~J{Q'\Q')^,n{Q'\Q') . (68) 

Eq. (|67|) is selfconsistent in the sense that the hamiltonian depends on the actual state 
I Q'{t) ). For example, if | Q'{t) ) is a single Slater determinant, but otherwise unrestricted, 
-^0 iQ'*j Q') is the Hartree-Fock hamiltonian ( Kerman and Koonin , 1976| ). Furthermore, the 
approximate Eq. is usually not a linear equation like the exact Schrodinger equation 
and therefore violates the superposition principle of quantum mechanics. 

Both, selfconsistency and nonlinearity are also common to classical molecular dynamics 
where the classical Hamilton function H{ri{t), ■ ■ ■ , pi{t) , ■ ■ ■) depends on the actual physical 
state. The nonlinearity of H{ri{t), ■ ■ ■ ,pi{t), ■ ■ ■) plays an important role for statistical 
properties like equilibration and chaotic behaviour, see Sec. 0. In quantum mechanics these 
properties can only be investigated by breaking the superposition principle by coarse graining 
or phase averaging. 



2. Quantum branching 

The time-dependent variational principle creates an approximate quantum dynamics in 
a manifold of trial states. Due to the severe restrictions it is very likely that areas of the 
Hilbert space are not reachable by the approximate equations of motion which would be 
visited by the solution of the Schrodinger equation. The idea for an improved model is to 
allow branching from one trajectory | Qi{t) ) to another one | Q'j(t) ), i.e. to jump in the 
parameter manifold with a certain probability. 

If the local deviation between the exact solution and the solution of the time- dependent 
variational principle, see Sec. |1I A 2| , is not so large, the operator 

= {H - q^^) 

= H-H^, (69) 
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which appears in Eq. (^) for the error, may be regarded as a perturbation, which causes 
transitions between the trial states | Qi{t) ), while each | Q'i{t) ) follows the approximate time 
evolution as given in Eqs. ( |65D and (|6^) with its selfconsistent hamiltonian Hq{Q'* ,Q'j). 
As in usual pertubation theory the transition \Q'i{t)) \Q'j{'^)) should be related to 
the amplitude ( Q'jit) \ AH \ Q[{t) ) ( |'l'ully| , |199(J| ; [I'opaler et a/.| , |1997| ; [Lacroix, (Jhomaz, and] 



Ayik| , [19^98, 1999|) . One would, however, like to work with probabilities instead of amplitudes. 



This should be possible if the physical system is in an energy regime with high level density 
where statistical arguments may be used. 



In literature one often employs phenomenological arguments (see Sec. ^) for random 
jumps in parameter space. The above sketched procedure which refers to AH has several 
advantages. The most important is selfconsistency, the transition depends on the actual 
situation and the choice of the trial states \Q'^{t)). If for example \Q'i{t)) is already a 
solution of the Schrodinger equation, AH \ Q'^it) ) = 0, no branching occures. 



Another important aspect is the breaking of symmetries ( |Colonna and Chomazj , |1998D 



Whereas the exact quantum state can retain dynamically conserved symmetries by linear 
superpositions of states, the nonlinearity and selfconsistency of the approximate hamiltonian 
(|55|) usually does not permit that. Take for example the mirror symmetry of a molecule or 
a nucleus which breaks into pieces. Each measured final channel which corresponds to one 
I ) does not possess mirror symmetry but the ensemble of all measured channels does. 
If the initial trial state \Q'{t = 0) ) has this symmetry and H does not break it explicitly 
the symmetry will be kept in | Q'(t) ). Here, quantum branching into a nonsymmetric trial 
state I Q'j ) and its mirrored counterpart | —Q'j ) with equal probabilities would resolve the 
problem of spurious cross channel coupling observed in selfconsistent approximations ([Griftin 



\et al\ , |198CI| ). The transitions may also allow quantum fluctuations to configurations which 
are close in Hilbert space but cannot be reached by the approximate time evolution, for 
example tunneling. 

Another noteworthy example of quantum branching is displayed in Fig. ^ which in a 
pictorial way shows four gaussian wave packets | qi{t) ) (see Eq. (^7|) ) forming a bound 
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FIG. 2. Sketch of a four particle cluster and a passing wave packet. Circles indicate 
half-density contours. The molecular dynamics equations of motion for the centroids 
and the widths of the packets can only influence the wave packet as a whole, either it 
passes by (b) or it is captured (c). The exact solution allows to split the packet into a 
captured and proceeding part (a) (see text). 

cluster and a wave packet (large circle) of an unbound particle which passes by. The time- 
dependent width of this wave packet has spread in coordinate space because it was moving 
freely for some time. Such situations occurr frequently in simulations. In the exact case 
a piece of this wave packet would stick to the four particle cluster indicating that, with 
a certain probability, a five particle cluster has been created. The remaining part of the 
wave packet would move on. Because the gaussian single-particle state does not permit this 
freedom the wave packet can only be bound to the cluster in total or escape. However, in 
order to become bound its width has to shrink to the size of the cluster otherwise the matrix 
element of the interaction with the other wave packets is too weak. It is obvious that beyond 
a certain size of the passing packet the overlap with the cluster will be too small to create a 
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sufficiently large generalized force on its width parameter to let it shrink. Therefore, small 
capture probabilities cannot be described. 

A branching procedure would in the above given example divide the time-evolution into 
a branch of a five particle cluster and one of a four particle cluster and a freely moving wave 
packet with the respective probabilities. In general one can say whenever the interaction 
varies strongly across the region of the wave packet additional shape degrees of freedom 
should be considered. Otherwise local actions are washed out and only mean field properties 
survive. 

But quantum branching faces also problems. In going from amplitudes to probabilities 
one should not violate conservation laws. For instance, one should not alter the energy 
distribution of the system. The mean value and the variance of the energy of the ensemble 
of trial states, which are populated during the branching should be the same as in the initial 
state. The same is true for other conserved quantities like total momentum or angular mo- 
mentum. It is quite conceivable that variational principles with the appropriate constraints 
could be helpful ( jBalian and Veneroni] , |1988| ; [Lacroix, Chomaz, and AyT^ , [1998, 19^ ). 



3. Approximation or new dynamical model? 

At this point some philosophical thoughts on the meaning of the results from the time- 
dependent variational principle seem to be appropriate. One could of course always claim 
that any restriction of the degrees of freedom in the trial state | Q' ) leads to approximations 
of the Schrodinger equation, the more restraint the trial state the worse the approximation. 



But the simple examples in Sec. |IIB| and those for antisymmetric states in the following 
Sec. m demand a more differentiated point of view. 

Let us exemplify this by restricting step by step the degrees of freedom in the general 
quantum state for a number of atoms. First we make a restricted ansatz for the trial state 
which contains only the coordinates and spins of nucleons and electrons and we disregard 
completely the internal quark and gluon degrees of freedom for the nucleons and the quanta of 
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the electromagnetic field. For low excitation energies this is certainly a good approximation. 

The next step is to neglect the internal degrees of freedom of the nuclei assuming that 
they are all in their ground states and to retain only their center of mass coordinates and the 
electron variables. But this is still too complex to solve the Schrodinger equation. Therefore 
we describe the cm. motion of the heavy nuclei by gaussian packets, which leads to the 
classical equations of motion (|53| ) for the nuclei with a Hamilton function 7i which couples to 
the quantal electrons. The state for the electrons may then be constrained to a single Slater 
determinant with no further restrictions on the single-particle states. Now we have arrived 
at the time-dependent Hartree Fock model which is a mean-field theory. The selfconsistent 
hamiltonian Hq is a one-body operator. 

But we can of course go further and disregard the internal degrees of freedom of the 
atoms or the molecules altogether and treat only their center of mass motion by means of 
fixed-width gaussians. At this point we obtain classical molecular dynamics with a Hamilton 
function 7i which contains two-body potentials between the molecules. 

If the molecules are big, may be junks of crystals, may be even stars hold together by 
gravity, their center of mass coordinates may be the adequate parameters for very narrow 
gaussians. The time-dependent variational principle will now provide Newton's equations 
for macroscopic objects. 

We don't believe that Newtons equations should be viewed as a bad approximation (bad 
because of the many constraints) of Quantum-Chromo-Dynamics where we started.^ They 
should also not be regarded as an approximation to time- dependent Hartree Fock, which 
formally one could plead for, because further constraints on the Slater determinant led there. 

The question is rather, which are the relevant degrees of freedom for a physical system. 
When those are identified the time-dependent variational principle provides a dynamical 



^The time-dependent variational principle can also be formulated in quantum field theory and it 
reduces to the nonrelativistic one discussed here. 
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model which is selfconsistent and has all properties of a Lagrange formalism, like existence 
of canonical variables, Noether's theorem, etc. 

Since the restriction of the trial state | Q' ) can be made in finer or coarser steps it is 
a matter of debate at which point one wants to speak about a new dynamical model and 
when of an approximation. 
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III. ANTISYMMETRIZATION 



The previous chapter demonstrated how to derive classical equations of motion from 
quantum mechanics by means of the time-dependent variational principle and localized 
single-particle states. It was also shown that there is no clear boundary between quantum 
and classical mechanics. This will be even more the case in this chapter where we are dealing 
with indistinguishable fermions. 

The natural generalization of the product state (^) (which led to classical molecular 
dynamics) to a trial state for identical fermions is the antisymmetrized product 



Q) = A I gi ) ® I g2 ) ® ■ ■ ■ ® I g^v ) (70) 
1 

iV! 



-L sgn(P) I gp(i) ) ® I gp(2) ) ® ■ ■ ■ ® | qp{^N) ) ■ 



all P 

The sum runs over all permutations P and sgn(P) is the sign of the permutation. The 
localized single-particle states | qi ) contain for spin-i particles also two-component spinors 



Xi) 



; X I ) = exp <{ — ) ®\xi) 1 bi = ri + iaiPi . (71) 



The spinor may for example be represented by the two complex parameters xj = ( T I X« ) 
and Xi = { l\xi)^ where 1 1 ) and | J, ) are the eigenstates of the spin operator Sz. 

The general Lagrangian C, as given in Eq. ([T^), calculated with the antisymmetrized trial 
state (^D|) provides through the Euler Lagrange equations (^T]) or ( P^ the desired molecular 
dynamics equations for fermions 

^Y.C,AQ*,Q)^. = ^^^^^^^■ (72) 

The hermitian matrix C^j^ {Q*,Q), Eq- (p3|), is not as simple as in the classical case, where 
C^iu = ^^lu, compare Eq. (^), but depends on all parameters contained in Q. 

Since the trial state is antisymmetric all consequences of the Pauli principle are incorpo- 
rated in the equations of motion (|72D. Cf^^, {Q*, Q), which is the second logarithmic derivative 
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of the determinant {Q\Q) = det{( Iqi)}, plays the role of a metric and will lead for ex- 
ample to large velocities q^, when the fermions get close to Pauli forbidden regions in phase 
space. In the energy Ti. {Q*, Q) = {Q\H \ Q)/{Q\Q) the Pauli principle causes exchange 
terms which induce for example additional momentum dependences. It should also be noted 
that the determinantal structure together with the nonorthogonality of the single-particle 
states result in an expression for the kinetic energy 

(Q\T\Q) JL 

7QfQ^=E(*lil*>0., (73) 

which contains a two-fold summation over all states (particles) compared to the single sum- 
mation in classical molecular dynamics, Eq. (^). Oik is the inverse of the overlap matrix 

0^')^^:={qk\qi); k,l = l,---,N. (74) 

The expectation value of a two-body operator like the interaction is a four-fold sum 

{Q\V\Q) 1 ( X 

{(Ikqi\v\q7nqn)[0nik0nl-0„u0nk)- (75) 



{Q\Q) 2 



k,Lm,n=l 



The generalized forces —-^'H iQ*^ Q) are therefore rather involved expressions which reflect 
the fact that the antisymmetrization ^ is a A^-body operation which correlates all particles 
simultaneously. 

The apparent large numerical effort led different authors to propose approximations 
which will be discussed in Sec. 0. 

First we explain in some length the two-body case because many of the new features 
concerning antisymmetrization and indistinguishability can be understood in this simple 
study. Furthermore, attempts to approximate the effects of the Pauli principle in molecular 
dynamics by so called Pauli potentials are based on considerations in two-body space. 

The consequences of antisymmetrization in many-body space, which are discussed in 



Sec. [IIIBI , are even more intricate. Depending on how much the single-particle states overlap 
the antisymmetrization can change the properties of the trial state completely. There is for 
example Fermi motion even if the gaussian single-particle states have no mean momentum. 
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A. Effects of antisymmetrization in two-body space 



1. Static considerations 

For two distinguishable particles the most simple trial state which leads in the proper 
limit to classical mechanics is the product of two gaussians \qi) ® | Q'2 ) • The corresponding 
state for indistinguishable fermions is the projection (ffOD onto the antisymmetric component 

I Q ) = A I gi ) ® I g2 ) = ^ { I gi ) ® I g2 ) - I g2 ) ® I gi )} . (76) 

In order to simplify the following discussion we do not use this Slater determinant but a 
trial state which separates center of mass and relative motion 

\Q)= \ qcm)® \ q) ■ (77) 
The center-of-mass wave function is parameterized by gem = B} as 

(l|gcm(t)) = exp |-^^^^||^| > m = m+zA{t)Pit) , (78) 

where X = ^ (xi + X2) is the center of mass coordinate and the parameter set qcm = {^7 B} 
contains the mean cm. position R the mean cm. momentum P combined in the complex 
parameter B and the complex width A. 

The wave packet \q) for the relative motion contains also the spins, thus g = 
{ a, b, Xi, X2}- Its parametrized form in coordinate space reads 

{x\qit))= exp|-^^^-^||xi(t))® |X2(t)) 
-ew{-^^^^}\X2{t))0\Xi{t)) 



(79) 



where small letters denote the relative coordinates and parameters. In relative coordinate 
space the exchange of the two particles, 1 ^ 2, is equivalent to the parity operation x = 
Xi —X2 ^ —X = X2 — Xi. 

The two-body wave packet (|77p is in general not a single Slater determinant because 
the center of mass motion for a Slater determinant does not separate if the single-particle 
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packets have different width parameters. In this section we choose the widths A{t) for the 
center of mass motion and a{t) for the relative packet to be independent in order to decouple 
the center-of-mass degree of freedom. It should also be noted that trial state ( |79D is a linear 
combination of total spin 5 = and S = 1. 

The relation between b and the more intuitive quantities "distance" r and "relative 
momentum" p is (time-dependence not indicated any longer) 



r — ajp and bf = a fip , 

a*b + ab , ^ lb- b* 
and p 



(80) 

a + a* I a + a* 

Here and throughout this section the real and imaginary part of complex numbers are 
denoted by the indices R and J, respectively. 

It is important to realize that r and p have their classical meaning only if the two particles 
are far away in phase space. We illustrate this effect of antisymmetrization in Fig. |^ where 
on the l.h.s. the relative wave function 

-3/4 ( 



'x\qd) 



2tt- 



a a 



exp 



X — r] 



+ ipx 



(82) 



a* + aj " [ 2 a 

is plotted along the r-direction. This trial state describes two distinguishable particles 
at relative distance r = {qd\x\qd). For a single gaussian the expectation value of x = 
x(l) — 2(2) is always equal to the parameter r independent of p and an. 

This is a typical example where Ehrenfest's theorem can be used to derive the classical 
equations of motion from quantum mechanics. The method is to replace the expectation 
values of the Heisenberg equation 



{qd\^i\qd) = {> 



2 fj,^ 



:) = -{<ld\k\qd) 



(83) 



for narrow wave packets by the mean values 

d ^ p 
dt /i 

Analogue for the relative momentum operator k 



(84) 



k{l) - k{2) 
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FIG. 3. Real (solid line) and imaginary part (dashed line) of the relative wave 
function for distinguishable particles (left column) and fermions (middle and right 
column) at different values of the parameters r = (r, 0,0), p = (p, 0, 0) and a = an. 
(d) measures the distance between particles, see Eq. (|88|). 
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.Qd\i 



V{x),k \qd) 



d 

ox 



(85) 



While ( |8^ ) is exact for | Q'd ), ( |86D is only an approximate expression which can be improved. 
However, the antisymmetrized packet 



{x\q) 



exp 



( X — r)^ 



+ ipx) — exp 



X + r) 



ipx 



(87) 



2a J \ 2a 

displayed for two values of r and p in the middle and right column leads to different results. 
First, for indistinguishable particles the operator x = 5(1) — 5(2) is not an observable 
anymore, because it is not symmetric with respect to particle exchange. Due to its negative 
parity the expectation value {(l\x\q)/{q\q) is always zero even for particles, which are far 
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away from each other. The same holds true for the relative momentum, {(i\k\q) / {q\q) = 
Therefore Ehrenfest's theorem is changed into a triviality because Eqs. ( PB|) and ( |53D become 
meaningless (0 = 0). The conclusion is that Ehrenfest's theorem cannot be used to derive 
classical equations of motion for indistinguishable particles. 

An observable which coincides with |r| at large distances is the root of the square radius 
minus the intrinsic width of the packet 

In the middle of Fig. ^ which is for p = 0, it becomes evident that the observable distance 
( d ) between two fermions with equal spin is larger than |r| and deviates most when <^ a^. 
The r.h.s. with \p\ = 2/ ^/or demonstrates that the distance ( d ) depends also on the relative 
momentum. For large relative momenta, 1, the Pauli principle is again less effective 

and (d) ^ \f\. The explicit expression for the mean square radius in terms of parameters 



IS 



k) .2 . 3|a|2 , lap 



{q\q} 2 an Ur 

where 

fex{^, Si2) = ZTT^ and C = = + qrp . 90) 

1 - e ^ 6i2 aR aR 

The first two terms in (|89|) are the same as for distinguishable particles. The remaining part 

is the exchange term which contains the quantity S, and the spin overlap = |( Xi I X2 )P- 

From Eq. (0) it is evident that may be regarded as a dimensionless measure for the 

distance in phase space where a/or is the length scale. For large ,^^4 the exchange vanishes 

and the Pauli principle is not active. For small ^ the distance between the fermions is always 

larger than |r| because the exchange term is positive definite. This effect is sometimes called 

"Pauli repulsion" , but this is not a force between the particles. 

The analogue expression for the distance in momentum space is 
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with the expectation value for the square momentum 



I I \ = P + 7^ + — /exU, ^12) 



(92) 



Again a term proportional to /ex(^, 'S'12) appears so that there is complete analogy between 
coordinate and momentum space. The exchange term is again positive and causes a "Pauli 
repulsion" in relative momentum. 

It is interesting to note that for a mininum uncertainty state | q ) with a/ = and equal 
spins (S'12 = 1) the sum of the observable distances in coordinate and momentum space 
fulfill 



{d) /aR+ {k) aji = ^' 



> 2 for a/ = . 



(93) 



Different from the measure the observable distance in phase space can never get smaller 



than 2, for details see subsection Sec. [Ill A4 . 

Expression (p^ tells also that the kinetic energy of relative motion 



T 



1 {(llffll) 



->2 



(94) 



2/i {q\q) 2yU ' AfiaR ' 2fiaR 

consists of the classical part p^/(2/i), the contribution from the uncertainty 3/(4/10^) and 
an additional potential fex{C,, S12) / {2fiaR) which depends on r,p, a and the spin overlap S'12. 
An analogue expression can be obtained for a spin-independent potential which is smooth 



{q\V(x) \q) 1 3|a|2 1 |aP 

V = Y ^ V{r) + i AV{r) 1^ + AV{r = 0) Su) 

{q\q) 6 2aR 6 an 



(95) 



This expression is exact if V{x) = Vq + V2 x^, and a good approximation if the Taylor 
expansion up to second order of V{x) around 5 = r is adequate within the range of the wave 
packet. Please note that due to rotational symmetry V{x) depends only on x^. 

Combining Eqs. (0) and (|95D yields a Hamilton function which splits into three parts 



n = T + V = U 



classical 



uncertainty 



+ v 



Pauli 



(96) 



+ 



4/iaij 4 ttR 
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Eq. (|9^) is a basis for the concept used by several authors in nuclear ([Wilets et aL\ , p.977 



IWilets et aT] , [T978| ; |Boal and GTosIi| , [T988| ; porso et al] , pHTt porso and Kandrupj , |I987| ; 



Peilert et al. 




1991 





\et al.\ , |1994| ; p<.lakow et al.\ , |1994| ; [Ebeling and Militzei| , |1997| ) to incorporate the uncertainty 



principle and the Pauli principle in order to extend classical molecular dynamics. A two- 
body potential J2i<j '^uncertainty ijij,Vij^(^) added to simulate the effects of the Heisenberg 
uncertainty principle and Y,i<j Vpauii(fij,Pij,a,Sij) is added to imitate the effects of the 
Pauli principle. The explicit form needs not to be the one given in (|96D, but it is usually 
adapted to the specific use. The method is quite successful in calculating energies ( Porso 
\al\ [19871 ; |Dorso and Randrup| , |1987| ), but we want to advise caution in using the Hamilton 



function (|9g) naively in equations of motion like 

^ -H and = 



d 



H . 



(97) 



dpj^ * dvi 

The reason is that and p.i are no longer canonical variables. As discussed earlier, although 
they still define the trial state | Q ) uniquely, they loose their intuitive meaning when the 
particles are indistinguishable. The operator x(i), "position of particle i" , is meaningless and 
( q I x(?) \q) / {q\q) is not r j ! The analogue holds for the momentum. Instead of postulating 
Eq. (|97|) we have to go back to the Lagrangian ( |T4|) and derive the equations of motion ([72|). 
It is evident that the matrix C will also be changed by the antisymmetrization. 



2. Center of mass motion 
Before discussing the relative motion in the antisymmetric case we consider first the 



center of mass motion. The Lagrange function (|I4D for the center of mass wave packet (|78| 
is given by 



r — r _'T 

'~'cm '~'0 cm cr 



(98) 



with 
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2 



(99) 



A* + A 



3 Ai 
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4 



B -B 

A* + A 



A* + A 



, \. \^ 3z /A A*^ 



P ■ : h total time derivative 



and the kinetic energy 



T 

r.r 



^ ^em I 2lf kcm ) 1 {B -Bf 



+ 



^2 
P 



(gem I gem) 2M + A 2M{A- + A) 2M AM Ar 

The Euler Lagrange equations yield 



5 = 0, A = — 
M 



(100) 



(101) 



or if one transforms B and B into R and P 

P 



P = 0, R 



M' M 



(102) 



In the center of mass wave function R{t) and P(t) have always the classical meaning of the 
mean center of mass position and momentum, respectively. Nevertheless, when the width 
parameter A is included as a dynamical variable, the wave packet | qcmit) ) = | A(t), B(t) ) 
is the exact solution of the Schrodinger equation. 



3. Relative motion 



The Lagrange function for the relative motion 



^q, q\ q, g*) = ^ ^^'7 , V'^^ - T ' " , , V (103) 
2 V {q\q) J {q\q) {q\q) 

features the antisymmetrization not only in the kinetic and potential energy as seen in 
Eqs. (p^ ) and (|95[) but also in the metric part Cq. 



For the sake of simplicity we treat in this subsection only equal, time-independent spins. 
In the general case the potential V{x^ 3(1), 3(2)) may of course depend on spin and the spin 
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degrees of freedom. The parameters ( t? i I Xi,2{t) ) would then also appear in Co,T and V 
and one would get equations of motion for them. For time-independent equal spins Cq is 
given by 



i (h 


* -h)(h +h) I 


/ 6 -6 \ 3 


2 


a* + a 4 


\a*+aj a* + a 



h-h -b-b 



[a —a) 



(104) 



b -b 



+ 



a* + a 



fexii, S12) + total time derivative 



^ J, 3 d/ 
p. r + h 

4aK 



p-r-r-p ^ {ag aj - aiaR-f -p cir) ^ o/ 



/ex (^5 S_ 



12) 



+ total time derivative . 

As expected, Co contains an exchange term, besides the terms for distinguishable particles, 
compare with Co cm in Eq- (|99|) - Although Co looks very complicated the Euler Lagrange 



equations (|72D 



iC 



dT 

d a* 



\ I 

+ 



av \ 

da* 



av 



(105) 



can be partially simplified because the free motion is, like for the cm. state or the case of 
distinguishable particles, again the exact solution of the Schrodinger equation 
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(106) 



and herewith 



dV \ 

da* 

dV 



(107) 



The complicated form of combines with the complicated derivatives ^) of the 

kinetic energy such that the simple result ( |106D is obtained. Although lengthy to calculate it 
is easy to understand. The reason is that both states, \a,b) and \a,—b), are exact solutions 
of the Schrodinger equation, they only differ in the initial conditions. The antisymmetric 
state \q) = \a,b) — \a,—b), defined in (|79D , is therefore also an exact solution. The 
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deeper reason is that both, the time derivative and the hamiltonian, commute with the 
antisymmetrization operator A in Eq. (|70D 



d 



d 



= z-A\^{t))-HA\<I^{t)) = A U- \^^f{t))-H\<^^{t)) 



dt 



dt 



(108) 



Therefore, if | ) is the exact solution of the Schrodinger equation, so is A | ). This 
general statement is not true if | ) is only an approximation. 

In the special case of a harmonic interaction, V{x) = Vq + ¥2^^, the contribution from 
the potential assumes a very simple form as well, because 



( d{x^) \ 
d a* 

V db* j 



y 2 a 6 y 



(109) 



The amazing result is that Eqs. ( |106| ) and ( |109| ) are the same for distinguishable particles 
and for indistinguishable particles, where is a complicated matrix depending on a*, a 
and h ,b. The exchange term in the expression for ( ) , Eq. (^), and for ( 5^ ) , Eq. (^9|), 
compensates for the different C. It is also interesting to see that for a harmonic interaction 
the parameters r and p obey the classical equations of motion 

A P 



fx 



p = -2 V2 r , 



(110) 



although the trial state | g ) is antisymmetric and describes two identical fermions. This 
result is, however, only obtained if the width a is at the same time a dynamical variable 
with the equation of motion 



/i 



i 2 1^2 a . 



;iii) 



The case, where a(t) = oq is supposed to be a positive time-independent number, leads to 
completely different results and is discussed in the following subsection. 



4- Relative motion with a time-independent width parameter 



This section investigates how the equations of motion change if the shape of the relative 
wave packet, Eq. (|79D , is restricted further by removing the width degree of freedom a{t) = 
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o_R(t) + iai{t) as a dynamical variable. For simplification only parallel spins are considered, 
i.e. 5*12 = 1. By setting aR{t) = Qq, aj{t) = and a{t) = in Eqs. (Q) and ( |104| ) we obtain 

Co =p-r + [p-r — r -pj — (112) 

r = ^ + ^ + ^^^, (113) 

2/i 4/iao 2/iao 1 — e ^ 

where 

^ = r^/ao+fao and C = Co-T-V. (114) 

The potential energy V is not given here explicitly and the spin dynamics is also not con- 
sidered for simplicity. The equations of motion for r and p are 

d dC dC 
= : — - or 

at dp op 

J. 2e"^ /a r{rr)/ao + p{pr)ao\ 2e~^ao 
-r - ^——^ ^——^ j + j^-^^ [ripp] - pirp)) 

= -|^(^'P) (115) 



and 



d dC dC 
= ■■ — or 

dt Of Or 



A , 2e « A r{rp)/ao+p{j)p)aQ\ 2e « / ^ ^ x 

^ + 1^ Y ) + (1 - e-^Ya, y^'^'^ - ""^''n 

= -^n{r,p). (116) 



For the reduced set of variables we are able to solve Eqs. ( |115| ) and ( 116 ) for r and p. The 
result is 

A f ^97^ ^dn\^ f r on ^dn\J 

= "i(0^ + "2(0 1 «o + p + — ^ - aoPis^ r (117) 



and 



ar / Vflo op or I 
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where and a2{^) are functions of = r^/ao +p'^clo and given by 



= (119) 



2e'^(l-e"g) ...n^ 

Since ai{^ ^ 1) = 1 and a2(^ 1) = one recognizes for ^ 1 Hamilton's equations 
of motion. ^ ^ 1 means that the two fermions are far from each other in phase space. In 
that hmit the identical fermions behave like classical distinguishable particles although their 
wave function is of course still antisymmetrized. When they get close in phase space, i.e. 
^<1, q;i(^)— s>^/2 and a2(0 3/^^ which means that the Hamilton-like parts in ( |117| ) 
and ( |118| ) vanish like ^/2 but the remaining parts increase like 3/^^. 



In this example one sees that for ^<;2 the equations of motion, which result from the 
parameterization (|7Sp with a{t) = a^, cannot be cast into Hamilton's form when r and p 
are regarded as canonical variables. To prove this statement let us suppose that a Hamilton 
function 'HpauuijiP) exists such that Eqs. (|117|) and ( |118| ) can be written as 

n = — and Pi = , 2 = 1, 2, 3 , (121) 

OPi OVi 

where i denotes the three spatial directions. Let us now disprove the existence of a function 
T^Pauii by calculating the mixed derivatives |^ and ^ which should add to zero if Eq. (|121|) 



is true. From the equations of motion ( |117|) and (|118| ) it is easy to verify that 



l^ + f^T^O. (122) 

OTk dpi 

This disproves the existence of a hamiltonian 'HpauU^fiP) as a function of r and p which 
would describe the fermionic dynamics derived from the ansatz (|7Sp with a{t) = oq for the 
wave function. 

It is, however, possible to find a pair of canonical variables (p, vf) 



P = \ l^^. r and n = J ]^^ p , (123) 
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which are non-hnear functions of the original variables ( ^araceno et al.\ , |1983D such 

that Cq in Eq. ( |112| ) assumes the canonical form 

£0 = ^ {n-i5-p-^) +\j^{rT>). (124) 



With these new variables the equations of motion aquire the form (|121D where TipauniPi 



T-C(f,p) is the total energy expressed in the new canonical variables. However, 'Hpauuij)^'^) 
cannot be expressed in a closed form because Eqs. ( p.23| ) cannot be solved for r(p, vf ) and 
p(p,7f). 

Equations (|121|) represent an approach often used in literature to incorporate the effects 
of the Pauli principle by means of a two-body interaction. Different forms of momentum- 
dependent potentials have been added to the classical hamiltonian or to the lagrangian 
( [Neumann and Ya\ |1994|) . Here we see that this method differs in several aspects from 
fermionic dynamics. Firstly, as discussed above the dynamical behaviour of r and p needs 
not be of Hamilton's type. Secondly, one should not replace (^(^)) simply by V{f) because 
even for narrow wave packets there is an exchange term which is not small for ^ < 1, see 
Eq. (^. 

The equations of motion simplify appreciably for a time-independent width parameter 
Oo, but the price to be paid is that free motion without interaction is not exact anymore. 
For < 1 the equations of motion differ essentially from the expected result r = and 
p = (see ( p.l7| ) and ( |118| )). The shape of the wave packet is restrained too much so that 



the particles scatter even if there is no interaction. 

There is however the appealing feature of the canonical variables (p, vf ) to exhibit a 
geometrically forbidden region in phase space. From (|123|) it is evident that p^/ao + vf^ao > 2. 
This is demonstrated in Fig. ^, where in a one-dimensional example several trajectories of the 
relative motion are shown for two freely moving fermions. The l.h.s. displays the trajectories, 
which are actually contours of constant kinetic energy, using r and p as dynamical variables, 
the r.h.s. shows the same trajectories but using p and tt. The empty area in the middle 
of the right-hand figure is the Pauli forbidden region of phase space. It corresponds to the 
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FIG. 4. Trajectories in phase space of one-dimensional free motion with fixed width. 
L.h.s.: variables r and p. R.h.s.: canonical variables p and tt. 

single point at the origin (r = and p = 0) in the left-hand figure. The trajectories which 
cross there, go around the circle in the canonical variables p and tt. 

Pauli potentials are usually chosen such that a pair of particles acquires a high energy in 
the forbidden region. One should however be aware that the kinetic energy at the boundary 
is finite, namely 5/(4/iao) in the three-dimensional case, see ( |113| ). The deficiency, that 



particles scatter even if there is no interaction, is also present in all Pauli potentials. In- 
clusion of the complex width a{t) as a dynamical variable cures this problem as has been 
demonstrated. 

At the end of this subsection the reader should not be left with the impression that 
the Pauli principle is a two-body effect. In fact antisymmetrization is a genuine A^-body 
correlation as will be discussed in the following. 



Resume 



1. The intuitive idea to include the Paul principle in a classical description by treating the 
exchange terms of the relative kinetic and potential energy of two identical fermions 
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as an additional "Pauli potential" which supplements the classical equations of motion 
cannot be supported. It is not correct to regard the parameters r and p as canonical 
variables if 

b-b* = \r + iap\'^ < 2\a* + a\ (125) 

and hence it is questionable whether the expectation value of the hamiltonian can be 
used as the Hamilton function in Hamilton's equations of motion. 

2. Heisenberg's quantum uncertainty refers to the variances ( ( 5^ ) — ( 2 ) ^) and {{k^) — 
( ^ ) ^) of the wave packet which are given by the width parameter a and are not related 
to or ]f. It is therefore open to doubt if inclusion of uncertainty into classical 
equations of motion can be achieved by a potential which depends on {ff, -Pij)- 
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B. Effects of antisymmetrization in many-body space 

1. Shell structure due to antisymmetrization 

It is not immediately obvious that an antisymmetrized product state like Eq. (|70|) includes 
shell-model features, like the nodal structure of single-particle orbits, because the states 
are localized in coordinate and momentum space. But due to the invariance of a Slater 
determinant under linearly independent transformations among the occupied single-particle 
states, after antisymmetrization, any set of single-particle states which is complete in the 
occupied phase space is as good as any other. This applies also to nonorthogonal states. 
To illustrate this we take four one-dimensional real gaussians with the same real width 
parameter Cq and zero mean momentum and displace them hj d = 0.75^/^0 (see l.h.s. of 
Fig. The one-body density can be written in terms of orthonormal states | ipm ) as 

N N 

p(i)= Oki {qi\ = E , (126) 

k,l=l m=l 

where the orthonormal eigenstates of p^^\ called natural orbits, are given by the following 
superposition of gaussians 

N 

l^m) = E kfc) (^^)fc- • (127) 
fe=l 



FIG. 5. Natural orbits. Antisymmetrization of four displaced gaussians (l.h.s.) 
leads to four occupied natural orbits (r.h.s.) which are almost harmonic oscillator 
states (dashed lines). 
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Oki is the inverse of the overlap matrix, {0~^)ki = {<ik\<ii)- They are displayed on the right 
hand side of Fig. ^ and compared to harmonic oscillator eigenstates (dashed lines). One 
observes that the occupied single-particle states | ipm ) consist of an s-, p- ,d- and an f-state, 
all very close to harmonic oscillator states. The difference between both sets can be made 
arbitrarily small by letting d/ y/oo approach zero. The f-state is for example essentially 
the first gaussian minus the second plus the third minus the fourth. All others are similar 
combinations. As already stressed in the two-body case in the previous section, when the 
wave packets overlap, ri and pi loose their classical meaning of position and momentum of 
particle /. In the limiting case d —>■ all ri —>■ and all — > and the harmonic oscilla- 
tor shells emerge. The distributions in coordinate and momentum space are the quantum 
mechanically correct ones of four spin polarized fermions in a harmonic oscillator. 

2. Fermi-Dirac distribution due to antisymmetrization 

A second example is illustrated in Fig. ^, where we consider 100 equally spaced gaussians 
in one dimension. Again all mean momenta are zero and the width oq is real. In the 
upper part of Fig. ^ the width y/d^ is 0.2 of the mean distance d so that the wave packets 
are well separated. Therefore the spatial density = {x\p^^^\x) and the momentum 
density pk = {k \ p^^^ \ k) are not changed by antisymmetrization. In the lower part the 
width has been increased to ^/oq = d. Without antisymmetrization (dash dotted line) the 
spatial density is uniform and the momentum distribution is that of a single packet. After 
antisymmetrization (full lines) one obtains the typical shell model oscillations in coordinate 
space and a Fermi distribution in momentum space. It is amazing to see how in Eq. (|126|) 
the superposition of gaussians by means of the inverse overlap matrices can create a fully 
occupied momentum state, see for example in Fig. |^ the lower right momentum distribution 
at /c = O.Skp, where the individual gaussians give practically zero probability to measure 
this momentum. We also calculated the eigenstates of the kinetic energy in the occupied 
space and got perfect sinusoidal waves. 
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FIG. 6. Densities in coordinate and momentum space. Full line with antisym- 
metrization, dash-dotted line without. Upper part: section of spatial density of hun- 
dred gaussians (not overlapping in coordinate space) and corresponding momentum 
distribution. Distributions with and without antisymmetrization are identical. Lower 
part: same as above but for overlapping gaussians. For details see text. 



These two examples illustrate nicely that Slater determinants of localized single-particle 
states with zero mean momentum can describe the harmonic oscillator shell model or even 
the Fermi motion of a gas of fermions in which plane waves are occupied up to the Fermi 
momentum. If one wants to simulate this effect by a "Pauli potential", disregarding the 
momentum distribution in each wave packet, the resulting ground state momentum distri- 
bution is unsatisfactory ( jDorso et al.[ |1987| ). The densities in coordinate and momentum 
space are just not given by the distributions of the ri and pi. 

At this point we would like to remind the reader of the different meaning of the three 
"momentum quantities" used in the article: 

1. momentum operator of particle / (not observable): k{l) = ^®---®A;®---Jl^, 
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2. momentum variable k as used in the momentum representation {k\q) and 

3. momentum parameter pi which characterizes the state \qi) = \ ri,pi, ai, xi )• 

Resume 

1. The coordinate distribution Px{x) is given by the observable 

= (QIQ) 

and not by the (eventually time averaged) distribution of the r^. Analogue the mo- 
mentum distribution is is given by 



{Q\ 






\Q) 


{Q\Q) 



Pkik) = (129) 



and not by the distribution of the Pi- All pi may be zero and there is still Fermi motion. 

2. Fermi motion is not random motion in and Pi- 

3. fi and Pi resume the classical meaning only when the system is so dilute in phase space 
that the Pauli principle has no consequences any more. 

3. Dynamical considerations 
The equations of motion for the most simple antisymmetric trial state (^) are, see 

Eq. 

^EC..(Q^Q)g. = U + |^. (130) 

Like in the two-body case, applying the inverse of C {Q*, Q) to the derivatives of the kinetic 
energy yields the simple result 



50 



ai 



hi 

0.2 
X2 

V • / 






i 

m 




V '■ J 



( 3V\ 

dV 

dy_ 

dV_ 

da* 

av 

all 

av 

V ''■ J 



:i311 



Although in the case of free motion ri and pi follow the classical path, this is once more the 
exact solution of the Schrodinger equation with a/ = i/m and bi and the spin parameters 
Xi being constant. The reason is of course that the many-body time evolution operator, 
which commutes with the antisymmetrization operator, is a product of single-particle time 
evolution operators and the free single-particle motion of gaussians is exact. 

The Pauli principle appears, like in the two-body case, in the interaction part in a twofold 
way. The interaction energy V {Q*, Q) has exchange terms and the metric {Q*, Q) cou- 
ples all generalized forces because it is in general not diagonal in the particle indices. 
Therefore the second term on the r.h.s. of Eq. (|131|) acts like a A^-body force. 

Looking at the structure of Eq. (|131|) it seems more natural to substitute the second 
term on the r.h.s. by a two-body interaction instead of considering the structure of the 
kinetic energy as done in Sec. [Ill A| . In nuclear physics Wilets and coworkers ([Wilets 
I], |1977| ; [Wilets et al\ |1978|) were the first who proposed a space and momentum dependent 



two-body "Pauli core" which was applied to atomic physics later ( [Wilets and Cohen[ , [19981 ). 
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IV. MODELS IN NUCLEAR AND ATOMIC PHYSICS 



In nuclear physics at non-relativistic energies the nucleons cannot be treated as classical 
particles on trajectories. Their phase space density is too large to ignore the Pauli principle. 
Therefore various molecular dynamics models which incorporate the Pauli principle on dif- 
ferent levels of sophistication are proposed in the literature. They will be discussed in more 
detail in the following subsections. 

Molecular dynamics calculations are also widespread in atomic physics to describe inter- 
acting atoms and molecules. But here the degrees of freedom are usually the classical cm. 
positions and velocities of the nuclei. The quantal electrons which move in the electric field 
of the atoms and provide the attraction are often treated by mean-field methods. There are 
also applications in which each electron is treated as an individual entity localized in phase 
space. Like for nucleons, if the density is high enough, they can become a degenerate Fermi 
gas for which a classical molecular dynamics picture cannot be applied. 



A. Antisymmetrized wave packets in nuclear physics 



The time-dependent Hartree Fock (TDHF) method ( pavies et al. , 1985 ) was expected to 



be suited for the description of colliding nuclei because stationary Hartree Fock calculations 
represent successfully ground state properties for all nuclei. The TDHF equations are ob- 
tained from the variational principle Eq. (§) when the single-particle states | q ) which form 
the single Slater determinant | Q ) are varied in an unrestricted way. However, it turned out 
that, even at low beam energies where the relative speed between the two colliding nuclei is 
small compared to the Fermi velocity and hence a mean field is always well established, the 
amount of fluctuations in the collective variables, like energy loss, scattering angle or mass 



1 


1985 




Balian 



and Veneroni| , |1992| pR^einhard and Suraud] , |1992| ; [Lacroix, Chomaz, and Ayik| , |1998, 1999D . 



Also the inclusion of collision terms, extended TDHF (see proceedings of "Time-dependend 
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Hartree-Fock and beyond" ( poeke and Reinhard] , |1982|) ), did not improve the situation. 
The reason is that the TDHF equations contain a common mean field which does not allow 
fiuctuations to grow. As discussed in Sec. [11 C 2| quantum branching into other Slater de- 
terminants with different mean fields is missing. Even the collision integral which induces 
fractional occupation of states and hence a mixing of Slater determinants does not change 
that because one common mean field is again calculated from this mixture. 

Calculations which treat the TDHF time-evolution for each member of a thermal en- 
semble with its own specific mean field ( [Knoll and Stracl^ , [1984[ ; [Knoll and Wu[ , [I988[) can 
describe large fiuctuations which develop during the expansion. The initial ensemble of 
Slater determinants, which could be envisaged as the result of quantum branching, was only 
assumed but not dynamically calculated. 

The shortcomings of a common mean field applies also to models which solve the Vlasov 
equation augmented by a collision integral like VUU (Vlasov Uehling Uhlenbeck) or BUU 
(Boltzmann Uehling Uhlenbeck), see for instance ( [Aichelin and Bertsch[ , |1985| ; [Stocker and 
[Greinei) , P^gB[ ; [Bertsch and Das Gupfa| , [Wolf et aL\ , [Gaitanos efoT] , jT^). 

Although no explicit antisymmetrization like in TDHF is present, Liouville's theorem and 
Pauli blocking in the collision term prevent an overoccupation of the single-particle phase 
space. The inherent lack of fiuctations in these one-body descriptions led to molecular 
dynamics models for fermions which will be discussed in the following subsections. 



1. Time- dependent a- cluster model 



The first molecular dynamics model which uses antisymmetrized many-body states of 
localized constituents in nuclear physics is an extension of the a-cluster model, which suc- 
cessfully characterizes a-particle nuclei, to the time-dependend case. In the model nuclei are 
represented as Slater determinants of wave packets for a-particles. The width parameter of 
the gaussian single-particle wave packets is chosen to be fixed ( [Caurier et al.[ [1982[ ; [Sara- 
ceno et oT] , [1983[ ) as well as time-dependent ( prozdz et aL\ , [19821 ). Even superpositions 
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of two Slater determinants of wave packets for a-particles with time-dependent width are 
investigated ( [Bauhotf et al\ , |1985| ). Since these models are applicable only for nuclei with 



pronounced a-substructure new models were developed in the late eighties, which address 
one wave packet to each single nucleon. 

2. Fermionic Molecular Dynamics - FMD 



The model of Fermionic Molecular Dynamics (FMD) was suggested in 1990 (|Feldmeie4 



19901 ; [Feldmeier et al\ , |1995| ; [Feldmeier and Schnack| , |1997|) in order to describe ground 



states of atomic nuclei and heavy-ion reactions in the energy regime below particle pro- 
duction. The many-body trial state of FMD is a Slater determinant \Q{t)) of single- 
particle gaussian wave packets | qi{t) ) where qi(t) denotes the set of single-particle param- 
eters qi{t) = {bi(t) , ai(t) , Xiit) , ^i} , bi(t) = ri(t) + iai{t)pi(t), which contains mean position, 
mean momentum and complex width. The spin degrees of freedom are represented by a 
spinor Ixiit)) (see Eq. ([TT])). The isospin part | ) is taken to be time-independent and 
identifies either a proton or a neutron. 

The equations of motion for all parameters are obtained from the variational principle as 
described in Sec. |l| (see Eq. (^)). The hamiltonian is an effective one, because the strong 
short range repulsion in the nucleon-nucleon interaction causes correlations which cannot 
be described by the trial state \Q{t)) ([Feldmeier et al\ , |1998| ). Also the spin correlations 
caused by the strong tensor force are only poorly represented in a single Slater determinant. 

In FMD the ground state is defined by the deepest minimum of the energy T-C{Q*, Q) = 
{Q \ H \ Q) . Since all generalized forces dTi/dgf = vanish in the minimum this state is 
stationary, all qi = 0. There is no "Fermi motion" in the parameters r;, pi, - ■ ■ . Ground state 
properties like binding energies and rms radii can be reproduced equally well with a variety 
of effective nucleon-nucleon interactions which differ mainly in their momentum dependence. 
But the intrinsic structure of the nuclei depends on the interaction. Superpositions of single- 
particle states as well as of Slater determinants can be used in order to obtain a more refined 
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description of nuclear structure, e.g. for halo nuclei ([Neff et al.[ |1999| ). 



FMD is able to model a variety of heavy-ion reactions ranging from fusion to dissipative 
reactions and multifragmentation ( [beldmeier and Schnack| , p,997| ). Different from TDHF 



fluctuations occur in these reactions but the results show also that initial correlations given 
by the intrinsic structure of the ground states play a major role in the simulation of frag- 
mentation reactions ( [Neff et al| , [1999|) . 

The time-dependent width parameters is an important non-classical degree of freedom 
( [Feldmeier et all |1995| ), especially to allow for evaporation of nucleons, a process that 
otherwise is strongly hindered, because each escaping wave packet takes away at least its 
zero-point energy. Inside a nucleus this zero-point energy is typically 10 MeV but evaporated 
nucleons have on the mean only 2 MeV kinetic energy. Therefore the packet has to spread 
during the evaporation process. 

Themo dynamic equilibrium properties can be determined in FMD by means of time 
averaging, see Sec. |V C| and Ref. ( [Schnack and I'eldmeierl , |1997| ) . 



As already explained in Sec. [11(J 2| the restricted parameterization leads to barriers which 
for the exact solution would not exist. Especially the splitting of wave packets is an im- 
portant source for quantum branching. The lack of this dynamical freedom is a serious 
hindrance of forming clusters ( [Kiderlen and Danielewic^ , [I996|) which is also observed in 
Antisymmetrized Molecular Dynamics and FMD investigations of spinodal decomposition 
of nuclear matter and quenches possibly important reaction mechanisms in multifragmen- 
tation reactions ( [Colonna and Choma2j , |199^ ). Further development in this direction is 



needed. 



3. Antisymmetrized Molecular Dynamics - AMD 

Antisymmetrized Molecular Dynamics (AMD) is similar to FMD with respect to the 
choice of the trial state but includes random branching between trial states. For details the 
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Horiuchi, p.996a| ; Pno and Horiuchi] , |1996b|) . AMD describes the nuclear many-body system 
by a Slater determinant \Q{t)) of gaussian wave packets chracterized by the parameter 



set qi(t) = {Zi(t),xi,C,i}, where in AMD notation the complex parameters Zi = —A=b, 



-^7=r; + i^J^pi denote the time-dependent centroids of the wave packets and xu the time- 
independent spin-isospin components which can be either proton or neutron, spin up or 
down. The width parameter = 2^;^ is real and time-independent and the same for all wave 
packets. 

The time evolution of the Zi(t) is determined by the time- dependent variational principle 
Eq. (^) which leads to the equations of motion given in Eq. (0). The Hamilton function 
Ti used in AMD is the expectation value of a hamiltonian H plus an additional term which 
removes the spurious zero-point energy of the center of mass wave packets for the different 
clusters. This cm. energy, which is of the order of 10 MeV for all clusters, is an artifact of 
all product states. 

The smooth variation of the Zi{t) due to the generalized forces dTi/dZl is supplemented 
by different stochastic forces. These can be regarded as a phenomenological ansatz for 
quantum branching between different trial states \Qj{t)) as discussed in Sec. [11 (J 2| . One 
branching procedure takes care of deviations caused by the short range repulsion between 
nucleons. For that a collision term is introduced which randomly changes the relative canon- 
ical momenta of a pair of wave packets. In order not to enter Pauli-forbidden regions in 
phase space approximate canonical variables are used which for two particles reduce to those 
discussed in Sec. [Ill A4| , compare Fig. ^. 

Another branching simulates the spreading and splitting of wave packets, which is an 
essential process for an adequate description of evaporation and absorption, but cannot be 
accomplished by the trial state ( phnishi and Randrup[ , [1995| ; [Ono and lloriucTii[ , [1996a[ ; 
[Ohnishi and Randrup| , [L997[ ; [Ohnishi and Randrupj , [1997] ). 

AMD is able to reproduce the essential properties of nuclear ground states (minima in 
H{Q*,Q)) and, when extended to trial states which use superpositions of single-particle 
states or Slater determinants, even the structure of halo nuclei, see e.g. ( [Kanada-En'yo el 
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al, |T995| ). 

Multifragmentation reactions are investigated for beam energies around the Fermi energy. 
Before comparing with experimental data from heavy-ion colhsions the resuh of a simulated 
collision is fed into a statistical decay program to account for long time processes ( pno| , 
Piihlhoferj 



Since the numerical effort of Antisymmetrized Molecular Dynamics as well as of FMD 
grows with A^^, for the calculation of systems with more than = 80 nucleons approxima- 
tions are needed. Recently the AMD group developed a "triple loop approximation" , which 
converts the fourfold sum of the potential energy (see Eq. (|75|)) into a threefold one, so that 
systems like Au + Au are feasible now (|Ono| , |1998| ). 
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B. Product states of wave packets — QMD 



Models which parameterize the many-fermion trial state by a simple product of gaus- 
sian wave packets are called Quantum Molecular Dynamics in nuclear physics. First QMD 
versions were invented in the eighties ( Aichelin and Stocker , 1986 ; Aichelin et al.[ 1987| ; Aiche 



[In , 1991 ; Khoa et~al. , 1992 ). They all employ a product state 



Qit) ) = I Qiit) ) ® I I2{t) ) ® ■ ■ ■ ® I g7v(t) ) (132) 



of single-particle states \qi{t)) = \fi{t),pi{t)) defined in Eq. p5|), where only the mean 
positions r;(t) and the mean momenta Pi(t) are time-dependent. The width is fixed and the 
same for all wave packets. 

The resulting equations of motion are the classical ones, given in Eq. (|53D . Also the 
interpretation of ri{t) and Pi{t) is purely classical and the particles are considered distin- 
guishable. This simplifies the collision term, which acts as a random force, and at higher 
energies the description of transitions from nucleons into resonances. 

All QMD versions use a collision term with Pauli blocking in addition to the classical 
dynamics. Some versions consider spin and isospin, others use nucleons with an average 
electric charge. Several QMD versions try to incorporate the Pauli principle by means of a 
Pauli potential that prevents nucleons of the same kind to come too close in phase space, see 
Sec. [Ill A| . Due to these simplifications QMD has the advantage that the numerical effort 
grows only with A^^ and thus allows the simulation of large systems. In addition all QMD 
versions use a statistical decay program for the long time dynamics. 

1. Versions 

QMD models are widely used in nuclear physics and exist in many versions. We would 
like to mention some of them and apologize for not mentioning all others. A more elaborate 



overview on QMD models is provided in Ref. ( [Hartnack et al.[ |1998|) . 
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1. Using experience with VUU/BUU models ( [Aichelin and Bertsch] , |1985|) one of the first 
QMD versions ( [Aichelin and Stockei] , |1986| ; [Aichelin] , |1991| ) was suggested. It exploited 



the trial state ( |132| ) only insofar that the interaction gets an effective range due to the 
folding with the wave packets. In any other respect the model propagates point par- 
ticles on classical trajectories. The zero-point energy originating from localization is 
omitted. As initial states random distributions of mean coordinates and momenta are 
taken according to the experimental ground-state density profile and binding energy. 
This distribution is however not the ground state of the model hamiltonian but an 
unstable excited state. The model ground state is highly over-bound. This QMD ver- 
sion does not distinguish between protons and neutrons, all nucleons carry an average 
charge. 

2. The IQMD version treats the isospin explicitly but is the same in all other respects. 
The version was designed for the analysis of collective flow and pions ( [Hartnack ei\ 

^, Mm. 



3. Another branch of the QMD evolution uses a Pauli potential and takes the trial state 
( |132| ) more seriously ( [Peilert et al.\ , |1992| ; [Konopka et al.\ , p.995| ). Since nucleons are 
kept apart in phase space by the Pauli potential ([Peilert et all [1991[) the minimum of 
the hamiltonian determines the nuclear ground state. This is not a mere theoretical 
beauty, but very important if one wants to investigate the survival of initial ground- 
state correlations in the final products. An ensemble of random initial states is not 
able to answer such questions. 



4. In the same spirit the Japanese QMD ( [Maruyama et al\ , [1992| ; [Niita et "al| , [1995| ; [Chiba 
et al] , [1996[ ; [Maruyama et all [1996[ ) is constructed. Supplementary the width param- 



eter which is time-independent in all other versions is chosen to be time- dependent 
here (EQMD). Analogous to AMD the cm. zero-point energies are subtracted from 
the hamiltonian, which is difficult because of the changing width. Some observations 
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in EQMD are similar to those in FMD. Especially evaporation processes are described 
much better than with fixed width. Also for fusion reactions the time-dependent width 
plays a major role because the fusion cross sections are too small with fixed width. 



5. Another version of QMD was developed in Copenhagen and called NMD ( Pondorf~ei 
xl], [19951 ). 



ai 



6. Many attempts have been undertaken to extend the applicability of QMD towards 
higher (relativistic) energies. Models are for instance RQMD ( ^orge et al.[ |1989| ; 
Lehmann et [1995t ISotH, |1995|) and UrQMD (|Bass et al.\, [19981). 



2. Decoupling of center of mass and relative motion 

A prominent problem of many-body trial states expressed in terms of single-particle 
quantities, irrespective whether they are antisymmetrized or not, is the center of mass 
motion which does not separate from the relative motion. One attempt to solve the problem 
is the construction of a trial state, where the single-particle width parameters are replaced 
by a width matrix Aki{t) ( |Kiderlen and Danielewic"z| , [19961 ) 

{xi,...,xn\ Q{t) ) = exp <^ -Y.(^k -rk{t)) Aki{t) {xi -ri{t)) + i^Pk{t)xk > | x) , (133) 

y k,i fe J 

\x) is a normalized spin-isospin state. The dynamical freedom of the matrix elements 
Aki{t) allows this state to factorize into cm. and intrinsic degrees of freedom for subgroups 
of particles. For two particles the advantage of ansatz ( [133| ) can be seen immediately 

(5i,52|Q) =exp|-Ae™(X-i?)%zPx| exp{-v4,eK2-^)^ + «p2} Ix)- (134) 

With the usual definitions of relative (5, r,p) and cm. {X,R,P) coordinates the width 
matrix-elements are related by Acm = 2(v4ii + A12), A^ei = {An — Ai2)/2, and An = A22- 
Independent of the relative motion, which can be in a bound state with width Arei, the 
cm. wave packet width Acm{t) can spread according to free motion. For product states the 
variance in the cm. coordinate is always connected to that of the relative motion. For the 
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product of two identical gaussian wave packets the relation is Acm{t) = 4:Arei{t) ( [Kiderlen 
|and Danielewicz| , |1996|) . 

The separation of internal and cm. variables should enhance fragment production which 
is otherwise suppressed due to the localization energy of the cm. motion. The authors could 
demonstrate a significant improvement in the description of light fragments. Unfortunately 
the proposed ansatz could not be extended to antisymmetrized states up to now. The trial 
state (|133D is also not flexible enough to describe the splitting of wave packets neccessary to 



model particle capture as discussed in connection with quantum branching in Sec. [11 C 2 



3. Approximate canonical variables - Pauli potential 



The historically first models that tried to describe fragmentation reactions on the basis 
of single-particle motion were classical models, see for instance ( [Bodmer and Pano4 |1977| ; 
Bodmer et al\ , |1980D . Problems arising from the fact, that classical particles neither obey 
the Heisenberg uncertainty relation nor the Pauli principle of identical fermions are attacked 
by introducing two-body interactions J2i<j '^uncertainty(^ij iPij ■, Ci) ( |Wilets et al.\ , |1977| ; |Wilets 
\et alj |1978| ) and J2i<j '^Pauuij'ij^ij^ ^t; Sij) ([Wilets et al\ |1977| ; [Wilets et a/.| , |1978| ; [Dorso ei 
olj |19871 ; porso and Randrupl, |1987| ; |Boal and Glosli] , |1988| ; [Peilert et a/.| , [I991| ; [Maruyama 
et a/.| , |1992| ; [Niita et~al\ , |1995| ; [Ebeling and Militzei) , |1997| ) which imitate the two quantum 



effects, see Eq. ( p6|) . 

The method is quite successful in calculating energies and reasonable single-particle 
occupations in momentum space for the free Fermi gas at finite temperatures ( Porso e\ 



I], |1987| ; Porso and Randrup] , |1987|) , but as already mentioned in Sec. pi A| , care should be 



taken in using these variables in Hamilton's equation of motion as if they were canonical. In 
addition classical models cannot correctly describe quantum statistical properties in general, 
like occupation numbers, mean energy, specific heat etc. 

Nevertheless, the idea of simulating Heisenberg uncertainty and Pauli exclusion principle 
by means of two-body interactions is still being used in nowadays apphcations of classical 
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dynamics to many-fermion problems, e.g. ( [Latora et al.\ , [1994^ [Wilets and Cohen] , |1998| ). 
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C. Atomic physics 



In atomic physics molecular dynamics applications are widespread because the de Broglie 
wave length of the atoms (molecules) is often much shorter than variations in the inter- 
molecular potential. Therefore a trial state for the center of mass coordinates, which is a 
product of well localized gaussians as discussed in Sec. |11 B 2|, is well suited 



Q'atom ) = I 5i ) ® • • ■ ® \Bn) with Bi = Ri + iAoPl 



(135) 



The electrons are much lighter so that quantum effects like Pauli principle and uncertainty 
are important. Their trial state might be thought of as an antisymmetrized many-body 
state I Q'^f, Q'atom. ) which depends on electronic degrees of freedom summarized in the set 
Q'ei = {q'o, Q'l, ?2, ■ ■ ■}, e.g. characterizing different orbits. It also depends on the variables 
of the atoms {Ri,Pi} which mark for example the phase space centers of the orbits. The 
dependence on Pi is usually neglected because the velocities of the electrons are much larger 
than those of the atoms. The total trial state is the product 



Q')= \ Q atom 



(136) 



Different from nuclear physics the hamiltonian is known. It can be well approximated by 
(spin and other relativistic effects neglected) 
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:i38) 



where capital and small letters denote atomic and electronic variables, respectively. From 
the Lagrange function Eq. (|^), one obtains with the appropriate approximations the well 
known Quantum Molecular Dynamics equations for atomic physics. The atomic variables Ri 
and Pi follow classical equations of motion (width Aq small) under the influence of electron 
potentials given by 
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+ ( Qel 'i Qatom \ el \ Qeh Qatom ) 



(139) 



The coupled equations of motion for the electronic degrees of freedom, Q'^i, are approximated 
at various levels of sophistication. 

These quantum molecular dynamics models are not subject of this review because they 
do not treat the electrons (fermions) as being localized in phase space. 



1. Product states of wave packets in atomic physics 



For simple Coulomb systems ( Juarez Barnes et al\ , |1993| ) and for hydrogen plasmas 
models are proposed which employ trial states of the type 



Q')= IQ'atom)® IQ'el) 



(140) 



where the electronic state | Q'^i ) consists of localized gaussian wave packets, Eq. (^5]). These 
models are QMD models in the sense of Sec. |1V B| , either with time-independent width and 
without Pauli potential ( [Ebeling and IVlilitzei) , |1997| ) or with Pauli potential ( [Ebeling and 
Schautzl , |19971 ) or with time-dependent width and Pauli potential (|Klakow et aLj |1994| ; 



Klakow et alj |1994|) . Only the electronic part is represented by a wave function, the protons 
of the hydrogen plasma move on classical trajectories. 

An interesting application, which is referred to in Sec. [VC2| , is the study of a plasma 
under extreme conditions, high temperature or pressure and phase transitions to the liquid 
and solid phase ( [Klakow et al\ , |1994| ; [Klakow et al.\ , [1994[ ). Other investigations focus on 
the degree of ionization of a partially ionised plasma (pbeling et al.[ [1996[ ; [Ebeling and 
Militzer[ , [T997[) . 
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2. Quantum branching 



When the Born-Oppenheimer approximation is vahd adiabatic energy surfaces can be 
calculated as 

^adiQatom^ Q'atom'^ ^el) = { ^eh Q'atom I X! 7^ ^7T7 + M^^l I Q'atom ) ? (141) 

1=1 i=i \Ri - X[l)\ 

where Vei = denotes the lowest energy state of the electrons under the influence of the 
charges Zi of static ions positioned at Ri. u^i = 1,2, ■ ■ ■ numerate the excited eigenstates of 
the electronic system, for example vibrational modes or particle-hole excitations. In principle 
one can set up an improved trail state as a linear combination 

\Q')=J2 Qu^^ I Q'atom,.^, ) ® I ^ei; Q'atom,u^, ) (142) 

and then try to solve the coupled equations which result from the variational principle for 
the complex amplitudes q^^i and the ion variables Qatom,u^i = 1-^; 'J = Ij ' ' ' i Uei = 

1, 2, ■ ■ ■}. From energy, momentum, and angular momentum conservation it is obvious that 
one has to have as many different trajectories {Ri P; for the ions as there are excited 
states which can be populated, because an inelastic excitation i/gi — >■ u'^i will change the ion 
trajectories accordingly. 

But these equations of motion are usually too complex to be solved numerically, one 
therefore introduces quantum branching (as discussed in Sec. 11 C 2| ) called quasiclassical 



trajectory surface hopping method. See for example references in ( p?opaler et al.\ , |1997D . 
Different approaches are tested against an accurate quantum dynamics calculation of a 
realistic system by Topaler et al. ([Topaler et al\ , |19971) . They find fair agreement between 



the quantum branching methods and results of the quantum equations which result from a 
trial state of type Eq. (|142|) . This is expected if electronic coherence is not important on 



time scales of the ion motion so that the hopping between the energy surfaces ( |141|) acts as 
a random Langevin force on the ion trajectories. 
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V. STATISTICAL PROPERTIES 



Molecular dynamics models are not only used to simulate nonequilibrium but also equilib- 
rium situations, especially when correlations require descriptions which go beyond mean-field 
approaches and quasi-particles. In the context of classical mechanics a vast literature exists 
(Hoover, 



, even for relativistic cases ( Belkacem et al. , 1998| ), in which equilibrium prop- 



erties are studied. So called thermostated time evolutions ( Nose , 1984 ; Hoover , |1985| ; Kus- 



nezov et~al. , 199U| ; Nose , 1991 ) , in which appropriate coupling to external degrees of freedom 



adjusts the temperature, are on firm grounds in classical mechanics and have been used 
successfully for equilibrium, e.g. to investigate classical spin systems, as well as for non- 
equilibrium, e.g. to study glass transitions. 

The classical procedures fail for quantum systems and especially for identical fermions 
when the phase-space density is not small anymore and the effects of the Pauli principle 
become important. No analogue to thermostated time evolutions exists yet for quantum 
systems. There are few attempts to infer thermodynamic properties from dynamical sim- 
ulations in quantum mechanics, some of them are given in ([Kusnezov| , |1993| ; [Ohnishi and 
Randrupl , |1993| ; [Klakow et al.\ , |1994| ; [Blaise et al.[ |1994| ; ^chnack and Feldmeiei) , |1996| ; Pno 
and Horiuchi] , |1996b| ; [Ebeling and Militzei) , |1997| ; [Ohnishi and Randrup] , |1997| ; phnishi and 
Randrupl , |1997| ; [Schnack and I'eldmeieil , |1997| [Schnack] , p.998| ). Usually one performs time- 
averages and relies on the ergodic assumption. The validity of this method depends of course 
crucially on the statistical and ergodic properties of the dynamical model. Following the 
mentioned articles one realizes that the matter is still under debate. 

Two general aspects are discussed in detail in the following sections. One concerns the 
thermostatic properties of a molecular dynamics model where the attribute thermostatic 
refers to the characteristics of the static canonical statistical operator, which determines the 
quantal partition function Z{T) = Tr(exp{— -^/T}). Once the trace is evaluated within a 
given model, its thermostatic properties can be deduced by standard methods like partial 
derivatives of In Z(T) with respect to temperature T or other parameters contained in the 
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Hamilton operator H. 

The other and even more important aspect, which is discussed in Sec. |V is the dynam- 
ical behaviour of a molecular dynamics model. For example a dissipative system which is 
initially far from equilibrium is expected to equilibrate towards the canonical ensemble. The 
simulation of such a system within the model provides a crucial test of its thermodynamic 
properties. Often one uses an ergodicity assumption, i.e. that time-averages are equivalent 
to ensemble averages, which should be verified. 



A. Thermostatics 

The question of theTmostatic properties can be condensed to the question, whether the 
set of model states \Q) is complete, i.e. able to span the many-particle Hilbert space. If 
that is the case the unit operator in A^-particle space can be written as 

(143) 

where d/i(Q) is a measure which depends on the parameter set Q because the trial states 
I Q ) are in general nonorthogonal. For fermions and gaussian wave packets the measure is 
derived in the following subsection. 

Once the completeness is shown the thermostatic relations have to be correct, provided 
the trace of the partion function 

Z(T) = Tr(exp{-^/T}) (144) 
{Q\exp{-H/T}\Q) 



^^^^^ (gig) 

is calculated with the trial states | g ) in quantum fashion. 

In the case of Fermionic Molecular Dynamics (FMD) and Antisymmetrized Molecular 
Dynamics (AMD) ( [Schnack and J^eldmeieij , [I^ ; pno et al\ , |T^^ ; |Ono et aL\ , |1992b| ) the 



antisymmetric many-body states form an overcomplete set and provide a full representation 
for the unit operator. As the calculation of the trace does not depend on the representation 
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all thermostatic properties like Fermi-Dirac distribution, specific heat, mean energy as a 
function of temperature etc. ought to be correct and fully quantal using FMD or AMD trial 
states. The difficult task is to calculate ( Q \ exp {—H^/T} | Q ), or ( Q | Bex^ {— H^/T} \ Q), 
where B is an observable. The trace integral can be evaluated by Monte-Carlo methods. 



1. Completeness relation with coherent states 

In the following it is shown that Slater determinants of coherent states span the whole 
Hilbert space for fermions. 

Coherent states 1 1 ) which are defined as the eigenstates of the harmonic oscillator 
destruction operator a = Wtt-x +i^/^k, 

1 / 3\ 

l\z) = z\z), /iHo = ^ (^a+a + -j (145) 
form an overcomplete set of states. They are the gaussian states defined in Eq. (HH) with 



^ ~ V 2a^^ ~'~ "^yi^P ~ Vv^^oo for a real width parameter qq. Their completeness relation 
reads in single-particle space 

^ r d^Rez d^lmz ipiai 
dV d^p I r, p)(r, p\ 



(27r)3 ( r, p I r, p ) 

where | z) labels the coherent states by their eigenvalue with respect to a^, and the phase 
space notation \r, p) labels the states by their expectations values of coordinate and 



momentum operators. Coherent states are extensively discussed in (|Klauder and Skager- 
itaml , [19851 ). 



Since we are dealing with fermions the spin degree of freedom has to be considered and 
consequently the resolution of unity changes to 

1(1)^ /dVd^ y MM (147) 
~ J (27r)3 4r {q\q) ' ^ ^ 

where the sum runs over the two magnetic quantum numbers m = ±i which are included 

in the set of parameters denoted by q (see Eq. (|7ll) ). 

68 



Proceeding one step further the unit operator in the antisymmetric part of the two- 
particle Hilbert space is the antisymmetric product of two single-particle unit operators 

1(2) = (iW^lW) (148) 



= (i'"®!'") i(l-£r2 

which may be expressed with antisymmetric two-particle states | gi, q2 )a as 

(2) ^ r dVi d^pi ^ r dVa d^P2 \qi, q2)a a{qi, 92 I 

J (27r)3 (27r)3 {qi\qi){q2\q2) ^ ^ 

where | gi, g2 )a := ^ ( ki ) ® 1^2 ) - 1^2)® | gi )) • 

Following this line the resolution of unity in the antisymmetric part of the A^-particle Hilbert 
space can be written as the projection of the A^-particle unit operator onto the antisymmetric 
subspace. Be \ Q) the unnormalized Slater determinant ([70|) of single-particle states | q ) 



IQ) = ]^Esgn(P)(|gp(i))®---® |gpw)) , (150) 

then the unit operator is 
with a measure 

mQ)-iQ\Qm ^^^,Jl^r Y,, (152) 

that accounts for antisymmetrization. In a sampling where the values of and are chosen 
according to Monte-Carlo methods this measure determines the probability to find the state 
I Q) belonging to the configuration Q = {ri,P]^,mi;r2,P2,m2; . . .} in Hilbert space. If for 
example two fermions with the same spin are close in r and p then this measure is small 
because the norm {Q\Q) = det{(gfc|g«)} will be small. Eq. ( |151| ) is very useful in calculating 
traces by means of Monte-Carlo sampling ( phnishi and Randrup| , |1993|) . 
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Coherent states are gaussian wave-packets with fixed width, but the real and imaginary 
part of the width, and a/, may also be integrated over appropriate ranges in order to get 
an improved coverage of the phase space when doing the sampling. 

|(i)_ [ dV d^p ^ f ^ f da/ \q){q\ , . 

J (27r)3 ^ k {q\q) ' ^ ' 

where and ^2/ denote the intervals the width a = aji + iaj is integrated over 

/ dan =Qr , / dai =Qi . 

Jur Jfii 

Since the width a in the completeness relation Eq. ( |146|) is arbitrary, the additional integra- 
tions in Eq. ( |153|) do not change the unit operator. 



In the case of nuclear physics one has two types of fermions, so that the measure for a 
system with neutrons and Z protons is 

Once the resolution of unity is given in terms of model states the partition function can 
be evaluated. 

It is interesting to note that the norm {Q\Q) in the measure cancels with the norm 
denominator in Eq. ( |144D . Provided the single-particle states are normalized, {qk\qk) = 1, 

and H does not depend on the spins, the partion function looks almost classical 
2N 
(2~) 

except that there is an operator in the exponent. Eq. (|155|) is still the exact quantum 
expression. There is no contradiction between the fact that the hamiltonian may have a 
discrete energy spectrum with gaps between the levels and the fact that the parameters 
and are integrated in a continous fashion ( [Schnack| , |1999| ). Only if the expectation value 



2^ r 

Z{T) = -—-^ j d\i d'p^ d\2 d'p2 ■ ■ ■ d'vN d^'pN ( g I exp {-H/T} \ Q) , (155) 



is moved up into the exponent things become wrong. The Pauli exclusion principle is fully 
taken care of by the antisymmetric state \Q)- If two particles with equal spin are at the 
same point in phase space | Q ) = and hence ( Q \ exp {—HjT} | Q ) = 0, so that forbidden 
states do not contribute to the partition function. 
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2. Example for many fermions 



In the following the above considerations are illustrated with the example of identical 
fermions in a common single-particle potential. Starting from the Hamilton operator 



TV 



1=1 



£(0 

2m 



+ ^(£(0) 



(156) 



the mean energy of the A^-fermion system can be derived from the partition function Z[T), 
Eq. ( |144| ), as the derivative with respect to T 

d 



JdfliQ) W(r) EZ,nOnm{T) [T'^{qm\ exp { - h/T} 



JdfliQ) W{T) 

where the two abbreviations VV(T) and 0'^{T) are introduced as 
W(T) 



(157) 



{Q 


exp{- 




Q) 


{Q 


Q) 



det[{qk\qi)^ 
{0-\T))^^ = {q,\ exp{-h/T} . 



:i58) 



(159) 



For free motion, i.e. v{x^) = 0, the matrix elements (|159|) are given by: 

' 27r alai 



exp 



2(4 + a; 



(160) 



.Ofe + «^ + ^ 

For the special case of the three-dimensional harmonic oscillator, i.e. f (x 
matrix elements can easily be calculated if one sets all = = l/{muj) ( [Schnack] , |1996| ): 



i mcj^x^, the 



\qk\ exp <; — — 



\qi) = ( — )' exp 
\muj/ 



muj r ^ 2 '^'^ ^* 



3iu 



2 bMl-e~-)-^\. (161) 



For other potentials the matrix elements ( |159D assume a more complicated form. 

Figure ^ shows the result of a Metropolis integration ([Metropolis et al.\ , |1953|) over and 
Pj^ of Eq. ( |157| ) for a system of four identical fermions with equal spins in a harmonic oscillator 
by open circles. The solid line displays the results when the traces are calculated by discrete 
sums over eigenstates of the many-body Hamilton operator ^ho ( pchmidt and Schnack] , 
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0.0 0.5 1.0 1.5 2.0 

T/mw) 

FIG. 7. Excitation energy as a function of temperature, calculated with antisym- 
metrized gaussian states (circles) or with eigenfunctions of the harmonic oscillator 
(solid line). The dotted line shows the classical result E* /N = 3T and the dashed line 
the result of quantum Boltzmann statistics. 



1998|) . For comparison the classical dependence is shown as a dotted line. The dashed line 
represents the result for the quantum Boltzmann case for distinguishable particles, were the 
trial state in Eq. ( |155| ) is a direct product of single-particle states. In both quantum cases 
the specific heat is vanishing at T = because of the finite energy gap between ground state 
and first excited state. 

As expected from the general argument that traces are independent of the representation, 
provided one uses a complete set, the numerical Metropolis integration with a sampling of 
10^ points in the 24-dimensional phase space of the continuous parameters rk,pf,, which 
specify the Slater determinants \Q), gives the same result as summing over eigenstates of 
the many-body hamiltonian. The continous integration is not in contradiction with Fermi 
statistics and the discrete spectrum of the hamiltonian, and the classically looking Eq. ( |155[ ) 



is fully quant al. 



Resume 
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1. The thermostatic properties of a model ought to be correct, if the set of model states 
I Q ) is complete, i.e. able to span the many-particle Hilbert space. 

2. In FMD and AMD the antisymmetric many-body states of single-particle gaussian 
wave packets form an overcomplete set and provide a full representation for the unit 
operator. 
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B. Thermodynamics 



In molecular dynamics the time evolution as given by the time-dependent variational 
principle, without a colhsion term or quantum branching, is deterministic. Given a state 
I Q{to) ) at an initial time to the system is described by the pure state | Q{t) ) at all earlier 
and later times. Therefore, like with the exact solution of the Schrodinger equation, thermal 
properties have to be obtained by coarse graining or time averaging. 

In this section time averaging is compared with the canonical statistical ensemble for a 
fermion system. If the system is ergodic both are equivalent and equilibrium properties can 
be evaluated by molecular dynamics simulations. For this the ergodic ensemble is defined 
by the statistical operator Rerg as 

te-ii) Jh {Q(t)\Q(t)) ^ ' 

Herewith the ergodic mean of an operator B is given by 

1 



( 5 ) I , « , - Tr (5=,. B) = ^lim j—, / d( 'TX^rA:::, ■ (163) 



tl 



{Qit) 


B 


\Q{t)) 


{Qit)\ 


Qit)) 



If the ergodic assumption is fulfilled, the statistical operator -Rcrg should depend only on 
( ^ ) , which is actually a constant of motion. Therefore the notation with the condition 
"(^)" in Eq. (|16D is used. 

Expectation values are well-defined with Eq. (|163D so that one can easily calculate exten- 



sive quantities like the excitation energy. But it is not obvious how an intensive thermody- 
namical quantity, such as the temperature, might be extracted from deterministic molecular 
dynamics with wave packets. 

In classical mechanics with momentum-independent interactions the partition function 

^classical(T) = J Y[ ^ ^(^2%)^^ ^classical , Pi , " ' ■)} (164) 

is a product of a term with the kinetic energy and a term containing the interactions. There- 
fore, a fit of the momentum distribution with a Boltzmann distribution, or the equipartition 

74 



theorem, can be used to determine the temperature T. For example, in simple formula- 
tions of the Nose-Hoover thermostat ( |Hoovei| , |1985| ; [Kusnezov et al\ , |199CI| ; |N ose| , |1991| ) the 
equipartition theorem serves as basic ingredient. 

In the quantum case Eq. ( |144| ) has to be employed which does not show this factorization. 
For an interacting finite system one cannot, in analogy to the Boltzmann case, fit a Fermi 
function to the momentum distribution to determine the temperature. An example is the 
ground state of a nucleus for which the momentum distribution has a smeared out Fermi 
edge due to the finite size and the two-body interaction and not because of temperature. 



1. Ergodic ensemble of fermions in a harmonic oscillator 

In this section the ideal gas of fermions in a common one-dimensional harmonic oscillator 
potential is used for demonstration. The hamiltonian ^ho is written in second quantization 
with the fermion creation operator (which creates a fermion in the n-th single-particle 
eigenstate of hno) as 

1=1 1=1 \ ^''^ ^ / n=0 

and the trial state \Q{t)) describes four fermions with equal spins. To test the fermionic 
nature of the dynamical evolution the ergodic ensemble averages of the occupation numbers, 
( c_,n£n ) \ ( H) ' ^'i' (HH) ' evaluated and compared with ( ( c^n^n of the canonical 
ensemble discussed in the previous section. 

The occupation numbers, which range from to 1 for fermions, are chosen on purpose 
to make clear from the beginning, that the equations of motion for the parameters Q{t) = 
{fi(t) , piit) , ■ ■ ■}, as given by Eq. (p2D, might be generalized Hamilton equations, but the 
observables always have to be calculated with the quantum state \Q{t)). If fk{t) and 
Pk{t) were taken as the particle coordinates in a classical picture, the question for the mean 
occupation number of the n-th eigenstate of the harmonic oscillator could even not be posed. 

As pointed out already, in Fermionic Molecular Dynamics (FMD) the time evolution of 
gaussian wave packets in a common oscillator is exact, and thus the occupation probabilities 
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of the eigenstates of the Hamilton operator do not change in time. In order to equihbrate 
the system a repulsive short-range two-body interaction is introduced. The strength of the 
interaction is chosen weak enough so that the ideal gas picture is approximately still valid. 

In the initial state | Q{to) ), which is far from equilibrium, three wave packets with a 
width of a(to) = ^/muj are put close to the origin at x(to) = {—d,0,d) with d = 0.5/-\/ma;, 
while the fourth packet is pulled away from the center. 

Figure || gives an impression of the chaotic time dependence of ( Q(t) \ £j£n | Q(t) ) for 
n = 0, 3 and 6. 




C0t/2JI 

FIG. 8. Occupation probabilities versus time - n = 0: circles, n = 3: squares, 
n = 6: triangles. 

The result of time averaging is seen in Fig. ^ (symbols) for four different initial 
displacements which correspond to four different excitation energies of the fermion sys- 
tem. To each case we assign a canonical ensemble which has the same mean energy, i.e. 
E* = {Huo - Eo)t ^ = {{Huo --Eo))|y The solid lines in Fig. |] show the cor- 
responding distributions of occupation probabilities for these canonical ensembles. Their 
temperatures T are also quoted in the figure. The one to one correspondence between the 
occupation probabilities of the ergodic ensemble and the ones of the canonical ensemble, 
which has the same mean energy ( ^ ) as the pure state, demonstrates that the system is 
ergodic and that the fermion molecular dynamics trajectory covers the many-body phase 
space according to Fermi-Dirac statistics. 
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This result is not trivial because, firstly, the system is very small, consisting of only four 
particles, and secondly, the equations of motion are approximated by FMD. 



2. Describing the system with fixed-width trial states 



As explained already in Sec. [V A| , both type of trial states with dynamical and fixed 
widths span the whole Hilbert space and thus the thermostatic properties are the same. 
Reducing the degrees of freedom by keeping all width parameters at a fixed value leads 



however to a different dynamical behavior as discussed in Sec. [Ill A| . The equations of motion 



are no longer exact solutions for the case of free motion. In a common oscillator the exact 
solution is obtained only if all ai = l/(mu;) because dai/dt is zero in that special case, see 
Sec. [IV A3|. If the width has a different value, spurious scattering occurs like in the case of 



free motion whenever two particles come too close in phase space. 

The left hand part of Fig. |10| displays for the very same system as in the previous 
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section the result of the time evolution but without the randomizing two-body interaction. 
If the widths are chosen to be a/ = l/{muj), the resulting exact time evolution is just a 
unitary transformation in the one-body space and the occupation probabilities are stationary 
(circles). But if the widths are taken as ai = 1.2/(mu;) the occupation probabilities change 
in time and the spurious scatterings equilibrate the system even without an interaction. The 
right hand part of Fig. |l^ shows the mean occupation probabilities when the interaction is 
switched on (triangles). Apparently the nature of the randomizing force is not relevant, it 
may even be a spurious force which originates from a too restricted trial state. 




FIG. 10. Occupation probabilities calculated in the ergodic ensemble using trial 
states with fixed widths (symbols) compared to those calculated in the canonical 
ensemble (solid line). L.h.s.: Without interaction, circles: ai = l/{muj), triangles: 
ai = 1.2/{mu>). R.h.s.: With interaction and ai = 1.2/(?nw). 



In Antisymmetrized Molecular Dynamics (AMD) ( |Ono et all |19924 |Ono et aL[ |1992b| ) 
trial states with time-independent widths are used, and as expected from the above simple 
example the thermodynamic properties of the model comply with Fermi-Dirac statistics. 
It would be interesting to see how collision term and branching influence the dynamical 
statistical properties of AMD. As the Pauli- blocking prescription is consistent with the 
AMD state we expect a faster equilibration due to the additional randomization. 
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3. Canonical and ergodic ensemble for distinguishable particles 

To complete the discussion in this section the fermions are replaced by distinguishable 
particles, i.e. the antisymmetrized many-body state is replaced by a product state of gaussian 
wave packets. In this case the relation between temperature T and excitation energy E* /N 
is known and given by E*/N = iu;[coth(^) — 1]. 



T 1 1 r 




5 10 15 20 

n 



FIG. 11. Mean occupation numbers for a product state (Boltzmann statistics). 
Ergodic ensemble (triangles), canonical ensemble (solid line), both for an excitation 
energy of E* /N = 2.05 MeV which corresponds to a temperature T = 5 MeV in the 
canonical ensemble. The dotted line shows the result for Fermi-Dirac statistics at the 
same temperature. 



The ergodic ensemble is again investigated at different energies and compared with the 
corresponding canonical ensembles with the same mean energies ( pchnack and Feldmeier , 



1996 ; Feldmeier and Schnack , 1997 ). The time evolution of the system exhibits also for this 



case ergodic behaviour at all excitation energies. As an example Fig. 11 is showing the case 
of E* /N = 2.05 MeV (which corresponds to T = 5 MeV) after a time averaging of about 
2000 periods. The ergodic ensemble (triangles) and the Boltzmann canonical ensemble (solid 
line) are the same within the size of the symbols. This means that the ergodic ensemble is 
again equivalent to the quantum canonical ensemble and not to the classical one, because 
for E* /N = 2.05 MeV and T = 5 MeV one is still in the quantal regime according to the 
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relation given above. The classical relation E*/N = T for the one-dimensional oscillator 
holds only for E*/N:^uj (here u = 8 MeV). 

However, since distinguishable particles are not affected by the Pauli principle, the oc- 
cupation numbers for the many-body ground state look quite different from those of the 
Fermi-Dirac distribution at the same temperature (dotted line in Fig. ^). 

Resume 

When discussing statistical properties of molecular dynamics with gaussian wave-packets 
(coherent states) one should always keep in mind that any observable or statistical weight 
has to be calculated with the trial state according to quantum mechanics. One should not 
fall into a completely classical approach, misled by the "classical" appearance of equations 
of motion or phase space integrals, which is due to the representation of the coherent states 
in terms of and Pj^. 

Statistical properties of molecular dynamics for fermions can be deduced from simulations 
of equilibrium situations but due to quantum effects a measure for intensive quantities like 
temperature is not readily available. 
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C. Thermal properties of interacting systems by time averaging 



In the previous section only noninteracting or weakly interacting systems are considered 
for which partition functions can often be calculated because the Hamilton operator is a 
one-body operator. If the interaction between the particles is strong enough, for example 
attractive enough to form self-bound many-body systems, the ideal gas picture is no longer 
valid and solid and liquid phases appear besides the vapor phase. Here the partition function 
Z{T) cannot be evaluated analytically anymore because it would amount to solving the full 
interacting many-body problem H \ ) = | \E'„ ) in the desired range of energies. 

Take for example a fermion system in a large spherical container with repulsive walls. At 
zero temperature the lowest eigenstate | \E'n=i ) describes a self-bound system in its internal 
ground state located at the center of the container. With increasing energies not only 
internal excitations and cm. motion of this drop occur but there is also the possibility to 
have two or more bound objects which are separated from each other and surrounded by 
vapor. They can be in different internal excitations with various cm. energies and vapor 
energies, all adding up to the total eigenenergy En- This means that the quantum number 
n enumerates not only the excited eigenstates but also the cm. degrees of freedom, the 
partition into different drops, and the fermion vapor state. The number of eigenstates in an 
energy interval increases rapidly with excitation energy. 

In principle one can deduce all thermodynamic properties from the level density, but it 
is obvious that for those complex and highly correlated states | ) it is impossible to solve 
the eigenvalue problem, neither analytically nor numerically. Therefore one tries to simulate 
correlated many-body systems by means of molecular dynamics (which describes the time 
evolution in an approximate way) and replaces the ensemble average by a time average. 
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1. FMD - the nuclear caloric curve 



As discussed in the section V B it can be shown that even a small system with only a 
few particles distributed over two different harmonic oscillator potentials is ergodic. So one 
can use one subsystem, for which the relation between excitation energy and temperature 
is known, as a thermometer to determine the temperature of the other. This idea has been 
used in FMD simulations of phase transitions in nuclei where the ideal Fermi gas picture does 
not apply because the nucleons are interacting by strong two-body forces. The construction 
of a thermometer is necessary because the temperature cannot be determined from the 
momentum distribution (pchnack and Feldmeiei) , |1997|) . In an interacting small fermion 
system the Fermi-edge is broadened not only due to temperature but also due to the finite 
size of the system and two-body correlations. 

A thermostat (large heat bath coupled to the system) is not advisable because phase 
transitions in small systems are recognized best in a micro-canonical situation where the 
energy distribution is within a narrow energy range and variations of the level density p{E), 
which indicate a phase transition, are not blurred by the Boltzmann factor e~^^'^ of the 
canonical ensemble. 

The thermometer consists of a quantum system of distinguishable particles which move 
in a common harmonic oscillator potential and interact with the nucleons. The nucleus itself 
is confined by a wide harmonic oscillator potential which serves as a containment. This is an 
important part of the setup because it keeps the evaporated nucleons (vapor) in the vicinity 
of the remaining liquid drop so that equilibration with the surrounding vapor can take place. 
The coupling between nucleons and thermometer particles is chosen to be weak, repulsive 
and of short range. It has to be as weak as possible in order not to influence the nuclear 
system too much. On the other hand it has to be strong enough to allow for reasonable 
equilibration times. 

The determination of the relation between excitation energy E* and temperature T 
(caloric curve) is done by time-averaging of the energy of the nucleonic system as well 



82 



as of the thermometer over a long period according to Eq. ( p.63|) . The time-averaged en- 
ergy of the thermometer Eth determines the temperature T through the known relation 
Eth/Nth = ^ujthcoth^ujth / 2T) for an ideal gas of distinguishable particles in a common har- 
monic oscillator potential with frequency Uth (quantum Boltzmann statistics). 
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FIG. 12. Caloric curve, l.h.s.: of ^^Mg, ^^Al and ^^^Ca at huj = 1 MeV, from 
( pchnack and Feldmeiei , 1997] ), r.h.s.: as determined by the Aladin group from the 
decay of spectator nuclei, from ( pochodzalla et 'al. , 1995| ). 



The resulting caloric curves for the nuclei ^^Mg, ^^Al and ^°Ca are displayed on the left 



hand side of Fig. [T^. All caloric curves clearly exhibit three different parts. Beginning at 
small excitation energies the temperature rises steeply with increasing energy as expected 
for a Fermi gas in the shell model. There the nucleons remain bound in the excited nucleus 
which behaves like a drop of liquid. At an excitation energy of 3 MeV per nucleon the curve 
flattens and stays almost constant up to about 11 MeV. This plateau at T ^ 5 to 6 MeV 
marks the coexistence phase where at low excitation energy one big drop is surrounded with 
low density vapor. With increasing energy the drop disolves more and more into vapor until 
all nucleons are unbound and the system has reached the vapor phase. The latent heat is 
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hence about 8 MeV at a pressure which is estimated to be close to zero. 



The caloric curve shown on the l.h.s. of Fig. has a striking similarity with the caloric 
curve determined experimentally from the fragmentation of colliding nuclei by the ALADIN 
group ( |Pochodzalla et al\ , |1995| ; Pr'ochodzalla| , |1997] ). Their results are displayed on the r.h.s. 



of the same figure. The position and the extension of the plateau agree quite well with the 
FMD calculation. Nevertheless, there are important differences between the experimental 
setup and the one used in the simulations. For further details see ( pochodzalla] , p.997| ; pchnack] 



and Feldmeiei) , |1997|) . 



2. Phase transitions of hydrogen plasma 

Hydrogen plasma under extreme conditions, high temperature or pressure is of great 
current interest, since it shows new structural, dynamical and electronic properties like 
orientational ordering, pressure induced metallization or changes in the vibronic spectra 
( IHemley and1Tao| , |T992t |Klakow et aJ] , |199|; |Knaup et a/] , [T999| ). 



One model employed in this context, is called wave-packet molecular dynamics. It is 
using gaussian packets with time-dependent widths and a Pauli potential, which is derived 
from the antisymmetrization of pairs of nucleons, see Sec. p^V C| and ( [Klakow et al.[ |1994D . 



256 protons and 256 electrons are distributed in a cubic box with periodic continuation in all 
directions. Since the protons are classical their temperature is simply given by their kinetic 
energy via the equipartition theorem. 

The equations of motion are followed over 6 ■ 10^^^ s. One observable which is sampled 
over this time period is the proton pair distribution function gpp{r) which is related to the 
probability of finding two particles at the distance r. The proton pair distribution reveals 



details about the binding and short range correlations in the system. In Fig. the pair 
distribution function gpp{r) is plotted for two temperatures. The peak at r = 1.3 a.u. signals 
that the protons are bound in H2 molecules. The authors find that due to medium effects 
this bond length is shifted to a smaller value compared to the free one which is r = 1.47 a.u. 
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FIG. 13. The proton pair distribution function gpp{r) at T = 300 K (top) and 
T = 3000 K (bottom). The soUd hne shows the simulation with 128 molecules, the 
dashed one with 32 molecules. Figure provided by P.-G. Reinhardt. 



in their model. At T = 300 K pronounced peaks can be seen around r = 4 a.u. which 
indicate a sohd structure of the H2 molecules. That these peaks are smoother at the higher 
temperature is interpreted as a signal that the system is now in the liquid phase. 

The authors also find that wave-packet molecular dynamics is very efficient for the dis- 
cussed problem and faster than the Car-Parrinello method ( Car and Parrinello , 1985 ) while 
providing comparable results ( Klakow et al. , 1994 ). 



85 



ACKNOWLEDGMENTS 



For supplying us with references or figures we would like to thank S. Bass, M. Belka- 
cem, Ph. Chomaz, P. Danielewicz, W. Ebeling, O. Knospe, Tomoyuki Maruyama, A. Ono, 
J. Randrup, E. Rass, P.-G. Reinhard, U. Saalmann, R. Schmidt and H. Stocker. We also 
thank R. Roth for carefully reading the manuscript. 



86 



REFERENCES 



Aichelin J., and G. Bertsch, 1985, "Numerical simulation of medium energy heavy ion reac- 
tions," Phys. Rev. C 53, 1730. 

Aichelin J., and H. Stocker, 1986, "Quantum Molecular Dynamics - a novel approach to 
iV-body correlations in heavy ion collisions," Phys. Lett. B 176, 14. 

Aichelin J., A. Rosenhauer, G. Peilert, H. Stocker, and W. Greiner, 1987, "Importance of 
momentum- dependent interactions for the extraction of the nuclear equation of state from 
high-energy heavy-ion collisions," Phys. Rev. Lett. 58, 1926. 

Aichelin J., 1991, "Quantum Molecular Dynamics - a dynamical microscopic A^-body ap- 
proach to investigate fragment formation and the nuclear equation of state in heavy ion 
coUisions," Phys. Rep. 202, 233. 

Arnol'd V.I., 1989, Mathematical Methods of Classical Mechanics, (Springer- Verlag, Berlin). 

Balian R., and M. Veneroni, 1988, "Static and dynamic variational principles for expectation 
values of observables," Ann. Phys. (N.Y.) 187, 29. 

Balian R., and M. Veneroni, 1992, "Correlations and fluctuations in static and dynamic 
mean-field approaches," Ann. Phys. (N.Y.) 216, 351. 

Bauhoff W., E. Caurier, B. Grammaticos, and M. Ploszajczak, 1985, "Description of light- 
ion collisions in the time-dependent cluster model," Phys. Rev. C 32, 1915. 

Bass S.A., M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, 
M. Hofmann, S. Hofmann, J. Konopka, G. Mao, L. Neise, S. Soff, C. Spieles, H. Weber, 
L.A. Winckelmann, H. Stocker, W. Greiner, Ch. Hartnack, J. Aichelin, and N. Amelin, 
1998, "Microscopic models for ultrarelativistic heavy ion collisions," Prog. Part. Nucl. 
Phys. 41, 225. 

Belkacem M., M. Brandstetter, S.A. Bass, M. Bleicher, L. Bravina, M.I. Gorenstein, 



87 



J. Konopka, L. Neise, C. Spieles, S. Soff, H. Weber, H. Stocker, and W. Greiner, 1998, 
" Equation of state, spectra and composition of hot and dense infinite liadronic matter in a 
microscopic transport model," Phys. Rev. C 58, 1727. 

Bertsch G.F., and S. Das Gupta, 1988, "A guide to microscopic models for intermediate 
energy heavy ion collisions," Phys. Rep. 160, 189. 

Blaise P., P. Durand, and O. Henri-Rousseau, 1994, "Irreversible evolution towards equi- 
librium of coupled quantum harmonic oscillators. A coarse grained approach," Physica A 
209, 51. 

Boal D.H., and J.N. Glosli, 1988, " Quasiparticle model for nuclear dynamics studies: Ground 
state properties," Phys. Rev. C 38, 1870. 

Bodmer A.R., and C.N. Panos, 1977, "Classical microscopic calculations of high-energy 
collisions of heavy ions," Phys. Rev. C 15, 1342. 

Bodmer A.R., C.N. Panos, and A.D. MacKellar, 1980, " Classical-equations-of-motion cal- 
culations of high-energy heavy-ion collisions," Phys. Rev. C 22, 1025. 

Bondorf J. P., D. Idler, and I.N. Mishustin, 1995, "Self-organization in expanding nuclear 
matter," Phys. Lett. B 359, 261. 

Broeckhove J., L. Lathouwers, E. Kesteloot, and P. van Leuven, 1988, "On the equivalence 
of the time- dependent variational principles," Chem. Phys. Lett. 149, 547. 

Broeckhove J., L. Lathouwers, and P. van Leuven, 1989, "Time-dependent variational prin- 
ciples and conservation laws in wavepacket dynamics," J. Phys. A: Math. Gen. 22, 4395. 

Car R., and M. Parrinello, 1985, "Unified approach for molecular dynamics and density 
functional theory," Phys. Rev. Lett. 55, 2471. 

Caurier E., B. Grammaticos, and T. Sami, 1982, "The time-dependent cluster model," Phys. 
Lett. B 109, 150. 



88 



Chiba S., O. Iwamoto, T. Fukahori, K. Niita, T. Maruyama, T. Maruyama, and A. Iwamoto, 
1996, "Analysis of proton-induced fragment production cross-sections by the quantum 
molecular dynamis plus statistical decay model," Phys. Rev. C 54, 285. 

Chomaz, Ph., M.Colonna, and A.Guarnera, 1996, "Spinodal deomposition of atomic nu- 
clei," in Advances in Nuclear Dynamics 2, edited by B. Arruada and G. Westfall, (Kluwer 
Academic / Plenum) . 

Colonna M., and Ph. Chomaz, 1998, "Spinodal decomposition in nuclear molecular dynam- 
ics," Phys. Lett. B 436, 1. 

Davies K.T.R., K.R.S. Devi, S.E. Koonin, and M.R. Strayer, 1985, "TDHF calculations of 
heavy-ion collisions," in Treatise of Heavy-Ion Science, Vol 3, edited by D.A. Bromley, 
(Plenum, New York), p. 3. 

Dorso C, S. Duarte, and J. Randrup, 1987, "Classical simulation of the Fermi gas," Phys. 
Lett. B 188, 287. 

Dorso C, and J. Randrup, 1987, "Classical simulation of nuclear systems," Phys. Lett. B 
215, 611. 

Dorso C, and J. Randrup, 1989, "Quasi-classical simulation of nuclear dynamics. Phase 
evolutuion of disassembling nuclei," Phys. Lett. B 232, 29. 

Drozdz S., J. Okolowicz, and M. Ploszajczak, 1982, "The time-dependent cluster theory - 
application to the a-a collision," Phys. Lett. B 109, 145. 

Drozdz S., M. Ploszajczak, and E. Caurier, 1986, "Variational approach to the Schrodinger 
dynamics in the Klauder's continuous representations," Ann. Phys. 171, 108. 

Ebeling W., A. Forster, and V.Yu. Podlipchuk, 1996, "Quantum wave-packets simulation of 
ionization processes in dense plasmas," Phys. Lett. A 218, 297. 

Ebeling W., and B. Militzer, 1997, "Quantum molecular dynamics of partially ionized plas- 



89 



mas," Phys. Lett. A 226, 298. 

Ebeling W., and F. Schautz, 1997, "Many particle simulations of the quantum electron gas 
using momentum-dependent potentials," Phys. Rev. E 56, 3498. 

Feldmeier H., 1990, "Fermionic molecular dynamics," Nucl. Phys. A 515, 147. 

Feldmeier H., K. Bieler, and J. Schnack, 1995, "Fermionic molecular dynamics for ground 
states and collisions of nuclei," Nucl. Phys. A 586, 493. 

Feldmeier H., and J. Schnack, 1997, "Fermionic molecular dynamics," Prog. Part. Nucl. 
Phys. 39, 393. 

Feldmeier H., T. Neff, R. Roth, and J. Schnack, 1998, "A unitary correlation operator 
method," Nucl. Phys. A 632, 61. 

J. Frenkel J., 1934, Wave mechanics, Advanced theory, (Claredon Press, Oxford), p 235. 

Gaitanos T., C. Fuchs, and H.H. Wolter, 1999, "Heavy ion collisions with nonequilibrium 
Dirac-Brueckner mean fields," Nucl. Phys. A 650, 97. 

Goeke K., and P.G. Reinhard, 1982, Eds., Time-dependend Hartree-Fock and beyond, Lecture 
Notes in Physics 171. 

Griffin J. J., P.G. Lichtner, and M. Dworzecka, 1980, " Time- dependent S'-matrix Hartree- 
Fock theory of complex reactions," Phys. Rev. G 21, 1351. 

Gerjuoy E., A.R.P. Rau, and Larry Spruch, 1983, "A unified formulation of the construction 
of variational principles," Rev. Mod. Phys. 55, 725. 

Hartnack G., L. Zhuxia, L. Neise, G. Peilert, A. Rosenhauer, H. Sorge, J. Aichelin, 
H. Stocker, and W. Greiner, 1989, "Quantum molecular dynamics: a microscopic model 
from UNILAG to GERN energies," Nucl. Phys. A 495, 303. 

Hartnack G., R.K. Puri, J. Aichelin, J. Konopka, S.A. Bass, H. Stocker, and W. Greiner, 



90 



1998, "Modelling the many-body dynamics of heavy ion collisions: Present status and future 
perspective," Eur. J. Phys. A 1, 151. 

Heller E.J., 1975, "Time-dependent approach to semiclassical dynamics," J. Chem. Phys. 
62, 1544. 

Hemley R.J., and H.K. Mao, 1992, "Anomalous low-frequency excitations in diamond-cell 
studies of hydrogen at megabar pressures," Phys. Lett. A 163, 429. 

W.G. Hoover, 1985, "Canonical dynamics: Equilibrium phase-space distributions," Phys. 
Rev. A 31, 1695. 

Hoover W.G., 1986, "Molecular dynamics," Lecture Notes in Physics 258, (Springer, Berlin) 

Kanada-En'yo Y., H. Horiuchi, and A. Ono, 1995, "Structure of Li and Be isotopes studied 
with antisymmetrized molecular dynamics," Phys. Rev. C 52, 628. 

Katz A., 1965, Classical Mechanics, Quantum Mechanics, Field Theory, (Academic Press, 
New York). 

Kerman A.K., and S.E. Koonin, 1976, "Hamiltonian formulation of time-dependent varia- 
tional principles for the many-body system," Ann. Phys. 100, 332. 

Khoa D.T., N. Ohtsuka, M.A. Matin, A. Faessler, S.W. Huang, E. Lehmann, and R.K. Puri, 
1992, "In-medium effects in the description of heavy-ion collisions with realistic NN inter- 
actions," Nucl. Phys. A 548, 102. 

Kiderlen D., and P. Danielewicz, 1996, "Fragments in gaussian wave-packet dynamics with 
and without correlations," Nucl. Phys. A 620, 346. 

Klakow D., C. Toepffer, and P.-G. Reinhard, 1994, "Hydrogen under extreme conditions," 
Phys. Lett. A 192, 55. 

Klakow D., C. Toepffer, and P.-G. Reinhard, 1994, "Semiclassical molecular dynamics for 
strongly coupled Coulomb systems," J. Chem. Phys. 101, 1. 



91 



Klauder J.R., and B.-S. Skagerstam, 1985, Coherent states, (World Scientific Publishing, 
Singapore) . 

Knaup, M., P.-G. Reinhard, and Ch. Toepffer, 1999, "Wave packet molecular dynamics 
simulations of hydrogen near the transition to a metallic fluid," Contrib. Plasma Phys. 39, 
57. 

Knoll J., and B. Strack, 1984, "The dynamics of the nuclear disassembly in a field-theoretical 
model at finite entropies," Phys. Lett. B 149, 45. 

Knoll J., and J. Wu, 1988, "Expansion dynamics and multifragmentation of a saturating 
system of fermions," Nucl. Phys. A 481, 173. 

Konopka J., H. Stocker, and W. Greiner, 1995, "On the impossibility of temperature ex- 
traction from heavy ion induced particle spectra," Nucl. Phys. A 583, 357c. 

Kramer P., and M. Saraceno, 1981, "Geometry of the time-dependent variational principle 
in quantum mechanics," Lecture Notes in Physics 140, (Springer, Berlin). 

Kusnezov D., A. Bulgac, and W. Bauer, 1990, "Canonical ensembles from chaos," Ann. of 
Phys. 204, 155. 

Kusnezov D., 1993, "Quantum ergodic wave functions from a thermal non-linear Schrodinger 
equation," Phys. Lett. A 184, 50. 

Lacroix D., Ph. Chomaz, and S. Ayik, 1998, 1999, "Quantal extension of mean-field dynam- 
ics," Proceedings XXXVI. Int. Winter Meeting on Nuclear Physics, Bormio/Italy p. 485, 
and "On the simulation of extended TDHF theory," Nucl. Phys. A 651, 369. 

Latora V., M. Belkacem, and A. Bonasera, 1994, "Dynamics of instabilities and intermit- 
tency," Phys. Rev. Lett. 73, 1765. 

Lehmann E., R.K. Puri, A. Faessler, G. Batko, and S.W. Huang, 1995, "Consequences of a 
covariant description of heavy ion reactions at intermediate energies," Phys. Rev. C 51, 



92 



2113. 

Maruyama Toshiki, A. Ohnishi, and H. Horiuchi, 1992, "Evolution of reaction mechanisms 
in the hght heavy-ion system," Phys. Rev. C 45, 2355. 

Maruyama Toshiki, A. Ono, A. Ohnishi, and H. Horiuchi, 1992, "Fragment mass distribution 
in intermediate energy heavy ion colhsions and the reaction time scale," Prog. Theor. Phys. 
87, 1367. 

Maruyama Toshiki, K. Niita, and A. Iwamoto, 1996, "Extension of quantum molecular 
dynamics and its application to heavy-ion collisions," Phys. Rev. C 53, 297. 

Metropolis M., A.W. Metropohs, and M.N. Rosenbluth, A.H. Teller, and E. Teller, 1953, 
"Equation of state calculations by fast computing machines," J. Chem. Phys. 21, 1087. 

Neff T., H. Feldmeier, R. Roth, and J. Schnack, 1999, "Realistic interactions and configura- 
tion mixing in fermionic molecular dynamics," Proc. Int. Workshop XXXVII in Hirschegg, 
ISSN 0720-8715, p. 283. 

Niita K., S. Chiba, T. Maruyama, T. Maruyama, H. Takada, T. Fukahori, Y. Nakahara, 
and A. Iwamoto, 1995, "Analysis of the (N,xN') reactions by quantum molecular dynamics 
plus statistical model," Phys. Rev. C 52, 2620. 

Neumann J.J. and G. Fai, 1994, "Classical lagrangian model of the Pauli principle," Phys. 
Lett. B 329, 419. 

Nose S., 1984, "A unified formulation of the constant temperature molecular dynamics 
methods," J. Chem. Phys. 81, 511. 

Nose S., 1991, "Constant temperature molecular dynamics methods," Prog, of Theor. Phys. 
Suppl. 103, 1. 

Ohnishi A., and J. Randrup, 1993, "Statistical properties of antisymmetrized molecular 
dynamics," Nucl. Phys. A 565, 474. 



93 



Ohnishi A., and J. Randrup, 1995, "Incorporation of quantum statistical features in molec- 
ular dynamics," Phys. Rev. Lett. 75, 596. 

Ohnishi A., and J. Randrup, 1997, "Inclusion of quantum fluctuations in wave packet dy- 
namics," Ann. Phys. 253, 279. 

Ohnishi A., and J. Randrup, 1997, "Quantum fluctuations affect the critical properties of 
noble gases," Phys. Rev. A 55, R3315. 

Ono A., H. Horiuchi, Toshiki Maruyama, and A. Onishi, 1992, "Fragment formation studied 
with antisymmetrized version of molecular dynamics with two nucleon collisions," Phys. 
Rev. Lett. 68, 2898. 

Ono A., H. Horiuchi, Toshiki Maruyama, and A. Onishi, 1992, "Antisymmetrized version of 
molecular dynamics with two-nucleon collisions and its application to heavy ion reactions," 
Prog. Theor. Phys. 87, 1185. 

Ono A., H. Horiuchi, Toshiki Maruyama, and A. Onishi, 1993, "Fragment formation studied 
with antisymmetrized version of molecular dynamics with two-nucleon collisions," Phys. 
Rev. C 47, 2652. 

Ono A., and H. Horiuchi, 1996, "Antisymmetrized molecular dynamics of wave packets with 
stochastic incorporation of Vlasov equation," Phys. Rev. C 53, 2958. 

Ono A., and H. Horiuchi, 1996, "Statistical properties of AMD for non-nucleon-emission 
and nucleon-emission process," Phys. Rev. C 53, 2341. 

Ono A., 1998, "Antisymmetrized molecular dynamics with quantum branching processes for 
collisions of heavy nuclei," Phys. Rev. C59, 853. 

Peilert G., J. Randrup, H. Stocker, and W. Greiner, 1991, "Clustering in nuclear matter at 
subsaturation densities," Phys. Lett. B 260, 271. 

Peilert G., J. Konopka, M. Blann, M.G. Mustafa, H. Stocker, and W. Greiner, 1992, "Dy- 



94 



namical treatment of Fermi motion in a microscopic description of heavy ion collisions," 
Phys. Rev. C 46, 1457. 

Pochodzalla J., T. Mohlenkamp, T. Rubehn, A. Scliiittauf, A. Worner, E. Zude, 
M. Begemann-Blaich, T. Blaich, C. Gross, H. Emling, A. Ferrero, G. Imme, I. lori, 
G.J. Kunde, W.D. Kunze, V. Lindenstruth, U. Lynen, A. Moroni, W.F.J. Miiller, B. Ocker, 
G. Raciti, H. Sann, C. Schwarz, W. Seidel, V. Serfling, J. Strotli, A. Trzcinski, W. Traut- 
mann, A. Tucholski, G. Verde, and B. Zwieglinski, 1995, "Probing the nuclear liquid-gas 
phase transition," Phys. Rev. Lett. 75, 1040. 

Pochodzalla J., 1997, "The search for the liquid-gas phase transition in nuclei," Prog. Part. 
Nucl. Phys. 39, 443. 

Piihlhofer F., 1977, "On the interpretation of evaporation residue mass distributions in 
heavy-ion induced fusion reactions," Nucl. Phys. A 280, 267. 

Reinhard P.-G., and E. Suraud, "Stochastic TDHF and the Boltzmann-Langevin equation," 
Ann. Phys. (NY) 216, 98. 

Saraceno M., P. Kramer, and F. Fernandez, 1983, "Time-dependent variational description 
of aa scattering," Nucl. Phys. A 405, 88. 

Schmidt H.-J., and J. Schnack, "Investigations on finite ideal quantum gases," Physica A 
260, 479. 

Schnack J., and H. Feldmeier, 1996, "Statistical properties of Fermionic Molecular Dynam- 
ics," Nucl. Phys. A 601, 181. 

Schnack J., and H. Feldmeier, 1997, "The nuclear liquid-gas phase transition within 
Fermionic Molecular Dynamics," Phys. Lett. B 409, 6. 

Schnack J., 1996, " Kurzreichweitige Korrelationen in der Fermionischen Molekulardynamik," 
dissertation. Technical University Darmstadt. 



95 



Schnack J., 1998, "Molecular dynamics investigations on a quantum system in a thermostat," 
Physica A 259, 49. 

Schnack J., 1999, "Thermodynamics of the harmonic oscillator using coherent states," Eu- 
rophys. Lett. 45, 647. 

Sorge H., H. Stocker, and W. Greiner, 1989, "Poincare invariant hamiltonian dynamics: 
Modelling multi-hadronic interactions in a phase space approach," Ann. Phys. (NY) 192, 
266. 

Sorge H., 1995, "Flavor production in Pb (160AGeV) on Pb collisions: Effect of color ropes 
and hadronic rescattering," Phys. Rev. C C52, 3291. 

Stocker H., and W. Greiner, 1986, "High-energy heavy ion collisions: probing the equation 
of state of highly excited hadronic matter," Phys. Rep. 137, 277. 

Suarez Barnes I.M., M. Nauenberg, M. Nockleby, and S. Tomsovic, 1993, " Semiclassical 
theory of quantum propagation: The Coulomb potential," Phys. Rev. Lett. 71, 1961. 

Topaler M.S., M.D. Hack, T.C. Allison, Yi-Ping Liu, S.L. Mielke, D.W. Schwenke, and 
D.G. Truhlar, 1997, "Validation of trajectory surface hopping methods against accurate 
quantum mechanical dynamics and semiclassical analysis of electronic-to- vibrational energy 
transfer," J. Chem. Phys. 106, 8699. 

Tsue Y., and Y. Fujiwara, 1991, " Time- dependent variational approach in terms of squeezed 
coherent states," Prog. Theor. Phys. 86, 443. 

TuUy J.C., 1990, "Molecular dynamics with electronic transitions," J. Chem. Phys. 93, 
1061. 

Wilets L., E.M. Henley, M. Kraft, and A.D. MacKellar, 1977, "Classical many-body model 
for heavy-ion collisions incorporating the Pauli principle," Nucl. Phys. A 282, 341. 

Wilets L., Y. Yariv, and R. Chestnut, 1978, "Classical many-body model for heavy-ion 



96 



collisions (II)," Nucl. Phys. A 301, 359. 

Wilets L., and J.S. Cohen, 1998, "Fermion molecular dynamics in atomic, molecular and 
optical physics," Contemporary Physics 39, 163. 

Wolf Gy., G. Batko, W. Gassing, U. Mosel, K. Niita, and M. Schafer, 1990, "Dilepton 
production in heavy-ion collisions," Nucl. Phys. A 517, 615. 



97 



